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F19628-76-C-0296.  Part  II,  dealing  with  a  3-D  ray  trace  study,  is  entitled:  "Ground 
Range  and  Bearing  Error  Determination  and  Display  for  an  OTH  Radar  Systems  with  an 
Arctic  Trough  Ionosphere". 

PM.  personnel  working  on  the  material  for  this  report  including  software 
development  include: 

T.B.  Barrett 
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The  report  is  incomplete  in  that  time  ran  out  before  theoretical  results  could 
be  properly  compared  with  experimental  results.  However,  we  believe  the  machinery 
exists,  largely  in  the  form  of  computer  software,  to  perform  the  calculations 
necessary  to  more  fully  correlate  theory  and  experiment. 
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Preface 


This  scientific  report  can  he  considered  to  be  mainly  on  the  "spectral  analysis 
of  radio  wave  amplitude  scintillations".  It  attempts  to  provide  both  a  pseudo- 
mathematical  scenario  (model)  as  the  basis  of  describing  and  "explaining"  certain 
observations,  and  what  might  be  best  characterized  as  technigues  for  relating  the 
model  with  observations.  There  are,  of  course,  many  physical  configurations  and 
parameters  which  give  rise  to  the  phenomenon  of  scintillation.  The  basics  are  a 
source  of  electromagnetic  radiation  (in  general,  light  waves  and  radio  waves  can  be 
studied  using  the  same  theory).  A  receiver-detector  (usually  assumed  in  the  far 
field  of  the  transmitter)  and  an  intervening  medium  which  is  at  least  partially  char¬ 
acterized  as  having  a  "random"  index  of  refraction.  No  one  model  (herein  referred 
to  as  general  scintillation  theory)  can,  from  a  practical  standpoint,  serve  to  des¬ 
cribe  all  realizations  of  the  above  physical  conf iguration  nor  does  it  make  any 
sense  to  try  to  model  the  general  situation.  This  report  attempts  to  deal  mainly 
with  what  has  become  a  "canonical  model"  (herein  referred  to  as  thin  screen  scintil¬ 
lation  theory),  i.e.,  transmitter  in  free  space  -  "thin"  perturbing  layer-receiver; 
both  transmitter  and  receiver  many  wavelengths  from  the  thin-layer.  Details  of  this 
model  and  its  range  of  applicability  to  reality  are  discussed  in  some  detail  in  this 
report. 

No  particular  attempt  is  made  in  this  report  to  justify  the  study  of  scintil¬ 
lation  in  this  particular  context.  It  is  merely  pointed  out  here  that  by  doing 
judicious  spectral  analysis  of  amplitude  scintillation  data,  it  is  possible,  assum¬ 
ing  the  validity  of  the  model,  to  partly  characterize  the  (remote)  scintillating 
layer.  From  this  standpoint,  the  analysis  of  radio  wave  scintillations  provide  a 
means  of  "probing"  certain  physical  characteristics  of  a  remote  medium  which  is  dif¬ 
ficult  and'expensive  to  study  by  more  direct  methods  (e.g.,  by  going  to  the  medium 
and  measuring  electron  density).  Spectral  analysis  of  scintillation  is  important  by 
itself,  however,  since  the  behavior  of  scintillation  "in  the  frequency  domain"  has 
important  consequencies  for  communications  systems  (though  this  aspect  of  scintilla¬ 
tion  analysis  is  not  discussed  further  here). 


The  study  of  general  scintillation  theory  requires  considerable  background  know 
ledge  in  diverse  areas  of  physics  and  math.  For  example,  electromagnetic  wave  propa 
gation  in  general  (including  optics  and  radio  waves  -  including  waves  in  a  complex 
ionosphere),  stochastic  process  theory,  the  theory  of  propagation  in  a  turbulent 
medium,  communications  theory,  partial  differential  equations,  and  for  results, 
various  areas  from  numerical  analysis  with  emphasis  on  computer  techniques. 

Any  report  reflects  first  of  all  the  author  "modus  operandi",  i.e.,  the  way  he 
likes  to  organize  and  think  about  things.  Secondly,  it  reflects  the  author's  under¬ 
standing  of  the  various  disciplines  which  together  provide  the  rational  for  the 
topic  under  study.  Finally  it  reflects  the  areas  of  "expertise"  where  the  author 
has  had  a  reasonable  amount  of  exposure  and  feels  comfortable.  This  is  reflected 
here  by  a  pre-disposition  towards  numerical  methods  and  examples. 

There  is,  a  vast  literature  on  general  scintillation  theory  and  on  what  is  basi 
cally  thin-screen  scintillation  theory.  Because  of  this,  the  subject  is  rife  with  a 
jargon  which  is  meaningful  mainly  to  those  who  are  scintillation  specialists  which 
the  author  is  not.  Thus,  terms  like  Fresnel  filtering,  S4  index,  etc.,  etc.,  are 
meaningful  to  scintillation  specialists  but  not  to  the  non-specialists.  Thus  this 
report  takes  some  care  to  define  and  justify  specialized  terms  which  this  author,  at 
least,  did  not  find  in  his  hoard  of  "basic"  knowledge.  Of  course,  what  is  basic  to 
one  person  is  new  to  another  so  there  is  a  limit  to  the  "didacity"  of  any  treatise 
or  report.  It  is  not  the  intent  of  this  report  to  be  a  survey  article  on  scintilla¬ 
tion  but  rather  to  provide  for  understanding  of  certain  aspects  of  scintillation 
under  conditions  that  are  made  clear  in  the  body  of  this  report. 

This  report  does  not  include  information  on  an  extensive  set  of  computer  soft¬ 
ware  which  was  developed  to  provide  for  data  reduction  (and  display).  Information 
on  this  software  is  contained  in  various  internal  PML  reports.  Since  it  is  particu¬ 
larly  germane  to  the  theme  of  this  Final  Report,  however,  an  extensive  discussion  of 
the  algorithms  used  for  experimental  spectral  analysis  is  included. 


There  is,  of  course,  no  standardization  of  symbols  for  scintillation  theory  as 
described  in  the  abundant  literature  on  the  subject.  For  a  single  report,  however, 
one  must  adhere  to  a  consistent  symbology.  We  have  chosen  to  base  ours  on  that 
which  appears  in  a  recent  "survey  book"  on  scintillation  theory  [Ishimaru  (1978)]. 

It  differs  slightly  mainly  in  the  use,  or  lack  of  them,  of  subscripts  and/or  super¬ 
scripts. 

e  (=dielectric  permitivity  or  just  plain  permitivity)  dielectric  constant 
which  describes  the  propagation  medium.  In  general  e  is  a  complex  tensor 
quantity  which  may  depend  on  field  strength.  In  this  report  it  is  always 
considered  to  be  a  real  scalar  quantity  independent  of  field  strength, 
n  wave  refractive  index  of  the  medium  (also  assumed  real)  unless  specified 
otherwise. 

The  relationship  e  =  n2  is  assumed  to  hold,  (see  note  1) 
a)  wave  radial  frequency 

A  wave  length  of  electromagnetic  field 

k  HL  =  wave  number 

A 

(subscript  o  denotes  wave  number  in  free  space) 
k  =  k0<n>  is  assumed  where  <n>  =  mean  wave  refractive  index. 

Note  that  k  is  always  assumed  to  be  a  non-fluctuating  quantity  -  it  is  a 
mean  wave  number  in  the  medium. 

U  a  (phasor)  component  of  the  electromagnetic  field  (U  is  generally  complex) 
ip  ln(U)  or  more  exactly  U  =  exp(tj>).  Since  U  is  complex,  \p  is. 

X,s  ip  =  X  +  is 

z  the  z-axis  is  usually  used  in  this  report  as  the  predominate  line  of  propa¬ 
gation  (typically  the  transmitter  and  receiver  are  placed  on  the  z-axis.) 
This  is  in  constrast  to  Ishimaru  where  he  uses  the  x-axis.  (Since  many 
current  radio  scintillation  problems  involve  looking  up  to  a  satellite, 
the  literature  on  radio  scintillation  commonly  uses  the  z-axis  while  light 
scintillation  most  commonly  uses  the  x-axis  since  the  propagation  of  laser 
beams  horizontally  through  an  atmosphere  is  the  typical  problem  at  light 
frequencies.) 


covariance  function  (=  <^(£  )  £{£2)>)  where  £(£)  is  usually  a  homogeneous 
random  process  and  £  is  a  point  in  Euclidian  space  with  dimension  2  or 
greater. 

Note  that  the  term  correlation  function  is  used  by  some  authors  while  we 

choose  to  reserve  the  term  for  <  £  -m,  £  -m>  where  m  is  the  process  mean. 

1  2 

If  the  process  has  zero  mean  (which  it  frequently  does),  there  is  no  dis¬ 
tinction  between  the  two  terms  and  they  are  used  interchangeably, 
spectral  density  function  associated  with  the  covariance  function  B.  It 
is  usually  assumed  that  $  and  B  are  Fourier  transforms  of  one  another. 

In  this  report,  the  stochastic  processes  are  considered  to  be  "locally" 

homogeneous.  That  is  to  say,  at  a  particular  point  in  space,  the  correla 

tion  function,  B(r  ,r  ),  associated  with  the  process  depends  only  on 
1  2 

r  -r  .  We  then  write  B(r  ,r  )  =  B(r)  =  B(x,y,z)  where  r  =  (x,y,z)  (local 
“1  “2  12 

coordinates).  If,  furthermore,  the  process  is  isotropic,  B  is  a  function 
of  |  r  -j^  |  which  is  symbolized  by  B(|£|),  j£|  =  +^+Z2  •  The  following 

relations  hold  between  the  correlation  function  B  and  spectral  density 
function  $  assuming  they  both  exist. 

B(£)  =  r  e1-*-  <t(K)  dK 

“  oo 

r°° 

=  J  cosK*£  <t  (K.)dK_  f or  real  random  fields 

since*(Kr_)  =  <j>(-K). 

*  e" B<l)dl  '(ip  C  cos<i'-)  B(^  - 

B ( | r  |)  =  B(r )  =  £jL  f  K$(K)sin(Kr)dK 
r  o 

♦dill)  =  *(K)  =  —■  /”  rB(r)sin(Kr)dr 

?TTZ  K  0 

For  homogeneous  and  isotropic  stochastic  fields,  there  is  an  associated 
(1-dimensional)  stochastic  process.  The  relationship  between  $  (K)  and 
the  spectral  density  function,  F(K),  associate  with  this  process  is: 


*(K)  =  -  1  iffil 

dK 

F(K)  =  2  tt  r  K${K)  dK 
o 

For  the  specific  case  of  "frozen-in"  turbulence,  the  following  relation¬ 
ship  exists  between  cj(K)  of  the  assumed  homogeneous  and  isotropic  stochas¬ 
tic  field  moving  with  mean  speed  v,  and  the  spectral  density  function  of 
the  "observed"  stationary  stochastic  process: 

S(M)  =  iL.  r  $(K)KdK 

V  |W|/v 

*<K)  =  Z±  S1  (Kv) 

2ir  K 

where 

Sl  (u)  =  d.Sfe») 

dw 

The  following  useful  relationships  esdst  between  B(jr),$(K)  and  B(x,pJ,  F(x,jc): 
B(x,e.)  =  /*0#F(x,jc)ei^-Cdic 

F(x.je)  =1  /  B(x,P)e"^'_p  dp 

Uir) 

F(x^c )  =  /°°  exp(iKxx)  $(K)dKx 
—  00 

4>(K)  =  JL  f°°  exp(-iKxx)  F(x,jc)dx 
2  TT  -°° 

If  B(x,P)  =  B(x,|g_|),  then  cylindrical  coordinates  can  be  used  to  obtain 
B(x . |jb |)  =  B ( x ,p )  =  2 it  /*  J0(kp)  ^ ( x » ic ) Kcitc 

F(x,|k|)  =  F(x,k)  =  A  /  JoO<P)  B(x,p)pdp 

2 IT  0 


£  a  generic  point  in  Euclidian  “geometric"  space  of  dimension  3. 

£  a  generic  point  in  Euclidian  “transform"  space  of  dimension  3. 

£  a  generic  point  in  Euclidian  "geometric"  space  of  dimension  2. 

£  a  generic  point  in  Euclidian  "transform"  space  of  dimension  2. 

re  classical  radius  of  the  electrons  =  ^*82  x  10*l3cm  (see  Jackson 

me2 

(1962),  p.  490).  This  constant  appears  frequently  in  the  literature  be¬ 
cause  of  the  relationship  between  the  index  of  refraction  n  and  the  elec¬ 
tron  density  N  in  a  plasma  medium  where  magnetic  field  and  collisions  can 
be  neglected.  This  relationship  is: 

n  =  (l-X)1/2  where  X  =  “P  where  u)D  is  the  plasma  frequency e  4ttN£2 
'  '  —  *  m 

(Gaussian  units)  and  ojis  the  radial  frequency  of  the  electromagnetic 
wave.  For  small  x  (n  close  to  unity)  we  obtain 


note  1  -  For  a  stochastic  propagation  media,  the  fluctuations  in  the  quan¬ 
tities  which  describe  the  media  are  usually  assumed  to  be  small.  If  x  is 
the  symbol  for  the  stochastic  quantity,  it  is  customary  to  write 
x  =  <x>(  1  +6x)  or  x  =  <x>(  1  +  x^  ) 

where  x  or  x  is  very  small  relative  to  1  and  x>  =  <Xj >  =  0.  The 
relationship  between  <5e  ~  8  n  follows  from 

£  =  <e>(l  +6e  ) 

n  =  <n>(  1  ♦  «n)  -  V-£-  -  gj 

Oo  £0 

If  <£>  *e0,  the  free  space  permitivity,  and  <n>  s  1,  we  find  1  +6n  = 

VT  +  5e  1  +  -  or  6  ~  26n.  Thus  it  is  permissible  to  equate,  for 

example,  second  moment  properties  of  the  fluctuating  component  of  permiti¬ 
vity  with  that  of  the  refractive  index,  except  for  a  factor  of  2. 

£  electric  field.  Specifically  £  =  £(£)  (£  is  a  complex  vector  quantity). 


6 
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X 

Y 

A 

<J> 
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The  received  (narrow  band)  signal  is  s(t)  basi¬ 
cally  R[U(t)eiut]  where  U(rJ  becomes  U(t)  due  to 
relative  motion  of  transmitter,  receiver  and  the 
gross  medium.  The  following  relations  hold: 


U ( t )  =  X(t)  +  i Y ( t )  =  A(t)e' 


s(t)  =  X(t)coscot  -  Y(t)sinojt 
A(t)  =  [X2(t)  +  Y2(t)]x/2 


<(>(t )  =  tan"1  Y(t)/X(t) 

a  quadratic  detector  produces  A2(t)  =  |U ( t ) | 2 ; 
a  linear  detector  produces  A(t)  =  |U(t)|  ; 

a  coherent  (linear)  detector  produces  X(t)  and  Y(t)  from  which  A ( t )  and/or 
<p  (t )  can  be  determined. 
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I nt roducti on 


In  this  report,  a  basis  for  explaining  certain  results  from  radio  wave 
amplitude  scintillation  measurements  are  reported  on.  The  measurements  were  made 
from  an  aircraft  flying  mainly  in  (magnetic)  equatorial  areas  with  the  transmitter 
located  in  a  satellite.  Some  of  the  results  of  these  experiments  have  been  reported 
on  elsewhere  (see  for  example  AFGL  (1978)).  These  reports  emphasize  some  of  the 
gross  features  of  the  observed  scintillation  and  their  relationship  to  other 
scintillation  measurements  from  ground  observations,  and  their  relationship  to  other 
observations  of  the  perturbing  medium.  The  emphasis  in  this  report  is  on  spectral 
analyses  of  the  measurements.  It  attempts  to  provide  a  model  in  which  to  imbed 
these  analyses  as  well  as  to  displaying  some  of  the  results  of  the  analyses. 

The  report  starts  out  by  showing  a  result  by  deWolf  (1975)  from  general  scintil¬ 
lation  theory  which  attempts  to  graphically  delimit  the  currently  accepted  models 
which  appear  to  be  useful  in  the  study  of  scintillation.  This  is  followed  by  a 
brief  description  of  the  various  theories  with  emphasis  on  the  thin-screen  theory 
since  this  latter  theory  is  the  one  chosen  to  be  used  here. 

The  "canonical"  thin-screen  scintillation  theory  is  then  described  for  the  gen¬ 
eral  configuration  and  relationships  between  the  spectra  which  characterize  the  per¬ 
turbing  medium  and  those  of  the  observed  amplitude  scintillations  are  presented. 
(These  results  are  mainly  those  developed  by  others,  e.g.,  Taylor  (1975).)  In 
general,  the  equations  which  result  are  not  amenable  to  solution  by  "analytical" 
means,  i.e.,  they  remain  as  complicated  integral  expressions  which  are  not  particu¬ 
larly  open  to  intuitive  interpretation.  However,  some  results  of  numerical  integra¬ 
tion  are  presented. 

The  canonical  model  is  then  "computerized",  so  to  speak,  for  the  specific  physi¬ 
cal  configurations  realized  by  the  satellite-thin  screen-ai rcraft  combination.  Here 
the  emphasis  is  on  providing  for  the  transition  between  the  predomi nantly  spatially 
varying  propagation  medium  (electron  density  fluctuations)  and  the  temporarily 
varying  radio  wave  amplitude  observations  (amplitude  scintillations). 


This  is  followed  by  a  brief  synopsis  of  the  art  of  spectral  analysis  with  empha¬ 
sis  on  the  method  used  for  results  presented  herein. 

In  summary,  this  report  provides  a  source  (or  reference  to  a  source)  of  the 
background  material  required  to  bridge  the  gap  between  the  theory  of  scinillation 
spectra  and  such  spectra  observed  experimentally.  In  particular,  the  report  leads, 
on  the  one  hand,  to  the  development  and  demonstration  of  a  computer  procedure  for 
producing  model  spectra  and,  on  the  other  hand,  to  the  development  of  data  reduction 
techniques  for  producing  the  experimental  spectra  to  be  compared  with  the  model 
ones. 


What  remains  to  be  done,  and  which  hopefully  will  be  documented  in  a  subsequent 
report,  is  to  attempt,  first  of  all,  to  find  good  fits  between  model  and  experiment 
and,  secondly,  to  determine  the  sensitivity  of  the  experimental  technique  to  changes 
in  various  model  parameters.  Thus,  typically,  one  might  investigate  how  (or  if) 
"tomograph  techniques"  (scan  the  scintillation  in  different  directions)  can  resolve 
various  layer  parameters  such  as  speed,  direction  of  motion,  anisotropy,  etc.  Going 
further,  how  do  the  multitude  of  parameters  tend  to  obfuscate  the  results?  Is  it 
really  possible  to  state  that  the  observed  spectrum  comes  from  a  "power  law  layer" 
as  opposed  to  a  "Gaussian  or  Kolmogorof  layer"? 
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Section  I  -  The  various  domains  of  scintillation  theory 


[From  D.A.  deWolf  (1975)] 


Before  summarizing  some  of  the  theoretical  methods  which  are  used  to  study 
general  scintillation  theory,  it  is  worthwhile  to  present  a  figure  which  attempts  to 
show  the  domains  of  appl icabi 1 ity  of  various  theories  (which  might  be  better  termed 
approximation  techniques).  In  addition,  this  figure  shows  approximately  the  differ¬ 
ence  between  the  "domains  of  interest"  of  most  optical  and  radio  wave  problems.  Per¬ 
haps  the  most  important  feature  of  this  figure  is  that  the  domains  are  2-dimension¬ 
al,  i.e.,  defined  by  two  scale  parameters.  One  parameter  (the  vertical  axis)  repre¬ 
sents  the  "number  of  mean  free  paths"  of  scattering  and  thus  delimits  "single  scat¬ 
ter"  regions  from  "multiple  scatter"  regions.  It  is,  of  course,  dependant  on  the 
statistical  properties  of  the  scattering  medium.  The  other  parameter  (a  length)  is 
merely  the  Fresnel  distance  associated  with  the  problem.  The  model  on  which  this 
diagram  is  based  is  very  simple,  i.e.,  a  plane  wave  perpendicularly  incident  on  a 
"slab  of  turbulence"  and  a  receiver  at  distance  L  within  the  slab.  The  turbulent 
slab  is  characterized  by  isotropic  and  homogeneous  turbulence  such  that  the  spectrum 
associated  with  refractive  index  variations  about  the  mean  level  is  the  Kolmogorov 
spectrum  given  by 


•  n(£)  =  *033  CfS) Kh1 1  2n  L^1  <|K|<  Zn  lo* 

=  0  |K|>  2u  lo1 

and  is  "undefined"  for  OK |K|<_  2-n  Lq1  •  Cn  is  the  structure  constant  which  is  related 
to  the  variance  of  the  refractive  index  variations  by  (see  Ishinaru  (1978)  p.  543) 


Cp  =  1.91  <6n^>  L^/j 
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The  mean  free  path  associated  with  this  spectrum  is 


a-J  =  [(. 033)(1.91)  <finO  \T0*h  /  K  K“11/J  dK]-1 

o 

where  K  =  |K  | 

Since  *n(K)  is  not  defined  between  0  and  21IL“o,  we  use  the  expendient  of 
integrating  between  21,!_o1  and  “  with  the  result 


I  t  K-b/J  dK  =  |(2")-VJ 
2"Lo 

Thus  a"1  ~  566  <<5170-'  k'*  Lq1 


The  now  classic  Kolmogorov  spectrum  for  refractive  index  fluctuations  caused  by 
turbulence  is  characterized  by  an  outer  scale  length  L0  and  an  inner  scale  length 
10,  the  idea  being  that  energy  is  transferred  within  the  "inertial  sub-range"  from 
large  turbulent  blobs  of  scale  size  L0  down  to  very  small  blobs  of  size  10  beyond 
which  energy  is  dissipated  by  molecular  collisions.  Such  a  model  seems  to  be  well 
confirmed  for  a  non-plasma  atmosphere  but  probably  is  not  valid  for  a  plasma  and  in 
particular,  for  the  ionosphere.  None  the  less,  the  essential  structure  of  the 
diagram  should  remain  valid. 


Figure  1.  Scintillation  approximation  regimes  (No  scale  but  imagine  logarithmic 

scales. ) 


U 


Note  that  many  of  the  approximations  used  in  scintillation  actually  were  devel¬ 
oped  for  scattering  problems  in  quantum  mechanics  since  the  Schroedinger  equation 
can  be  made  identical  to  the  scalar  v/ave  equation.  The  Moliere  approximation  is 
valid  for  small  angle  scattering  (see,  for  example,  Modesitt  (1971)  and  associated 
references). 

It  should  be  emphasized  that  this  diagram  applies  only  to  the  simple  model  as 
stated  and  the  thin  phase  screen  which  we  choose  to  emphasize  later  is  not  really  in 
this  diagram  (perhaps  other  dimensions  should  be  added  to  take  account  of  more  geome¬ 
tries).  None  the  less,  it  serves  as  an  introduction  to  some  of  the  approximations 
which  are  briefly  discussed  next. 


Section  II  -  Various  approximat ions  for  sci nti 1 1  at  ion  theory 


General  (wave)  scintillation  theory  is  imbedded  in  a  more  inclusive  theory 
called  wave  scattering  in  a  stochastic  medium.  This  theory  has  at  least  two  main 
branches  (transport  theory  might  be  considered  to  be  a  third):  discrete  scattering 
theory,  where  the  scatterers  usually  have  well  defined  boundaries  and  are  imbedded 
in  a  homogeneous  medium;  continuum  scattering  where  the  parameters  of  the  medium 
vary  continuously  but  stochastically  in  space  and/or  time.  An  overview  of  much  of 
wave  scattering  theory  is  presented  in  Ishimaru  (1978). 

The  continuum  theory  is  appropriate  in  the  study  of  line  of  sight  radio  wave 
propagation  through  a  layer  of  refractive  fluctuations.  This  theory  in  turn  has 
developed  along  several  paths.  Two  which  predominate  the  study  of  radio  wave  propa¬ 
gation  phenomena  are  the  "scattering  cross  section  approach"  and  the  "scalar  wave 
equation"  approach.  In  the  first  method,  a  stochastic  scattering  cross  section  is 
defined  over  a  volume  of  propagation  space  and  the  statistical  characteristics  of 
the  observed  field  (usually  far  from  the  scattering  volume)  are  obtained  directly 
from  formulations  of  the  scattering  cross  section.  This  method  is  a  generalization 
of  a  similar  procedure  used  for  discrete  scatterers.  (See,  for  example,  Ishimaru 
(1978),  Chapter  2  and  Chapter  16.) 

In  the  second  method,  which  is  generally  used  for  line  of  sight  propagation, 
the  starting  point  is  a  stochastic  differential  equation.  The  differential  operator 
may  operate  on  either  a  vector  or  scalar  field  but  because  of  obvious  computational 
considerations,  most  useful  results  which  have  been  obtained  are' for  scalar  fields. 
In  recent  years,  a  vast  literature  has  appeared  on  various  approximate  solutions  to 
the  stochastic  scalar  wave  equation.  This  literature  has  its  practical  origins  in 
the  area  of  laser  communications  systems  on  the  one  hand  and  satellite  radio  wave 
communication  systems  on  the  other  hand.  This  theory  is  generally  referred  to  as 
"wave  propagation  in  a  turbulent  medium"  wherein  it  is  usually  implied  that  propaga¬ 
tion  is  essentially  in  the  forward  direction  as  opposed  to  "scatter  propagation" 
wherein  the  communication  path  is  distinctly  bent.  (Gross  refractive  effects  may 
occur  in  "wave  propagation"  however.)  As  the  name  implies,  the  fluctuations  or 
"irregularities"  in  the  propagation  medium  are  due  primarily  to  turbulence.  Thus, 
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turbulence  theory  must  be  invoked  to  explain  the  statistical  characteristics  of  the 
fluctuations  (of  the  refractive  index,  for  example). 

The  remainder  of  this  section  is  devoted  mainly  to  surnnarizing  pertinent 
results  (or  approximations)  from  this  theory  -  referred  to  below  simply  as  stochas¬ 
tic  wave  theory.  Before  doing  so,  however,  it  is  pertinent  to  point  out  that  sto¬ 
chastic  wave  theory  (as  mentioned  previously)  is  a  subset  of  the  mathematical  theory 
of  stochastic  differentia?  and  integral  equations,  or  more  simply,  statistical 
continuum  theories.  In  Beran  (1968),  some  useful  generalities  concerning  this 
theory  are  presented  with  examples  of  solution  methods  in  different  physical 
regimes.  Beran  points  out  some  of  the  general  difficulties  in  solving  these  equa¬ 
tions.  Because  of  these  difficulties,  the  equations  for  the  stochastic  fields  them¬ 
selves  are  not  solved  but  rather  equations  for  various  statistical  moments  are  solv¬ 
ed  for  (typically  using  various  approximation  methods).  The  method  most  commonly 
used  in  stochastic  wave  theory  is  to  obtain  an  integral  equation  ''solution"  for  the 
wave  equation  and  to  then  determine  the  various  moments. 

The  starting  point  for  stochastic  wave  theory  is  the  stochastic  differential 
equation  for  electric  field  (in  general,  _E  is  a  complex  vector  quantity) 


E  +  ko  E  -  2V  (_Jn  •  E)  =  0  (2-1) 


—  _  n  ” 

See,  for  example,  Tatarski  (1961)  or  Ishimaru  (1978).  Particular  note  should  be 
made  of  the  fact  that  the  refractive  index,  n,  is  a  scalar  function  of  position  (and 
possibly  time)  only,  so  that  the  medium  is  assumed  to  be  isotropic  in  the  non-statis- 
tical  sense  of  the  term.  (The  fluctuations  of  refractive  index  may  be  anisotropic 
in  the  statistical  use  of  the  word,  i.e.,  the  function  Bfr^.rJ  =  <n(£A) ,n(£^)>  may 
be  a  function  of  r  s  |_ri  -  r *  |  only.)  Under  certain  conditions  (usually  assumed), 
the  last  term  in  equation  (2-1)  can  be  neglected.  Physically  this  means  that  polari¬ 
zation  effects  can  be  neglected  since  there  is  no  coupling  between  the  vector  charac¬ 
ter  of  the  electric  field  and  the  vector  electron  density  gradient.  The  simplifying 
assumption  is  that  A«  10  where  10  is  the  "scale  size"  of  the  stochastic  medium. 
Typically,  this  scale  size  is  defined  by  Bn(r)=  C  exp(-r/l0)  assuming  that  the 
correlation  function  for  n,  Bn,  has  the  above  form.  Other  forms  for  Bn  can  usually 
be  used  to  define  a  similar  scale  size.  With  this  approximation  (2-1)  can  be  simpli¬ 
fied  to  the  (scalar)  Helmholz  equation  (see  e.g..  Born  and  deWolf  (1964)  or  Jackson 
(1962)). 

(V*  +  ko  n^)U  =  0  (2-2) 

where  U  is  a  component  of  the  electric  field  vector  transverse  to  the  mean  propaga¬ 
tion  path.  A  "solution"  to  this  equation  (actually  the  integral  equation  form)  is 
found  using  the  outward  radiating  free  space  Green's  function  (e.g.,  see  Frisch 
(1968)). 

U  =  G<>(r,o)  -  k^  f  GO(r,r1)M(rl)U(ri)dV1  (2-3) 

V 

where 

k^M  =  ko(n^  -  <n>^)  (2-3a) 
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GOfr.r1) 


exp ( ik0  Jr-r1  j) 
ll'l1 1 


( 2- 3b ) 


In  order  to  truly  solve  such  an  expression  for  U  it  is  necessary  to  put  bounds 
on  the  region  of  integration  V  and  to  approximate  U  (the  unknown)  within  the  inte¬ 
gral.  Various  formal  solution  methods  to  2-3  are  given  in  the  above  reference.  For 
many  applications  the  Born  approximation,  sometimes  referred  to  as  the  "single- 
scattering"  approximation,  is  made  whereby  the  unknown  field  component,  U,  is  replac¬ 
ed  by  the  free  space  field  U0.  Another  technique  is  to  assume  the  form  U(r)  = 
W(r)eikz  for  a  wave  propagating  in  the  positive  z  direction  and  where  W(rJ  is  a 
function  which  varies,  at  most,  slowly  in  the  z-direction.  Then  the  Helmholz  equa¬ 
tion  is  written  as: 


2ikdW  +  V^w  +  k2  fin2W  =  0  (2-5) 

dz 


where 


kgn2  =  k2  (1+  fin2  ) 


Then,  assuming  the  above  mentioned  scaling  condition,  A«l0,  it  can  be  shown  that 


|k4W  |  »  | 

dz  dz2 


and  the  following  approximate  parabolic  equation  for  W  is  obtained: 


2ik^W  +  (4l  +  £_  )W  +  k2  fin2W  =  0  (2-6) 

dz  dx2  dy2 
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Another  approximation,  the  Rytov  approximation,  is  obtained  by  first  making  the 
substitution  U(r)  =  exp(v(£))  and  writing  y  =  where  y0  is  the  free  space 

component  of  v  and  Vj  is  that  due  to  the  irregularities.  The  equation  for  is: 


(v^  +  k^)U0  li-i  =  [v +  k^dn]U0 


where  U0  =  or  exp(v0(_r))  is  the  free  space  “component"  of  the  field  and 

dn  =  (1  +  <5n)z-l.  The  Rytov  approximation  ignores  terms  in  vy-j  with  the  result  that 
the  following  integral  equation  for  y-j  can  be  derived  (see,  for  example,  Ishimaru 
(1978)  or  Yeh  and  Liu  (1972)): 


Vi(r)=^i  f  (i°(£,£1)dn(ri)U0(ri)dV1  (2-7) 

U0(£)  v 


where  G°  is  defined  above. 

Equations  2-4  and  2-7  represent  "first  order"  solutions  for  the  scattered 
field.  They  are  "single  scatter"  solutions  in  that  within  the  scattering  medium, 
the  incident  field  is  always  the  free  space  field.  The  parabolic  equation  (2-6),  on 
the  other  hand,  can  be  used  for  "multiple  scatter"  or  "strong  scatter"  conditions. 

In  this  report  the  emphasis  is  on  the  parabolic  equation  approach  and,  in  particu¬ 
lar,  on  the  "thin  phase  screen"  approximation  to  the  parabolic  equation  to  be 
derived  later.  For  the  sake  of  completeness,  however,  some  results  from  the  Born 
and  Rytov  approximations  are  included. 
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Section  III  -  Second  moment  properties  of  the  observed  field 
for  various  approximations. 

The  emphasis  of  this  report  is  on  displaying  and  explaining  the  second  moment 
properties  of  fluctuations  of  the  observed  field  as  determined  from  scintillation 
ineasurements  (observations  of  the  "received  signal",  s  -  see  glossary).  These  sec¬ 
ond  moment  properties  can  equivalently  be  presented  in  terms  of  covariance  functions 
or  spectral  density  functions. 

For  the  Rytov  approximations,  we  have 


k/ 

*1  "  *  -  *o  a  ^ 


GO  dn  U0  dV1 


while  for  the  Born  approximation: 


Ui  =  U^o  =  ki  , 

U0  U0  U0  y 


GO  dn  U0  dV1 . 


Thus  the  second  moment  properties,  for  example,  of  ’•'•j  and 


Ui 

Tr~  are  identical  in  so 
uo 


far  as  the  assumed  approximations  apply.  As  is  shov/n  later  R i 3  is  equated  to  the 
signal  amplitude  in  db  (i.e.,  log  of  A).  Thus  consider  the  Rytov  approximation 


( 2-7). 


ev(r)  =  evo(r)  +  ^ i ( r )  =  u0(r)e^l  (£) 


vi(r)  =  X(r)  +  iS(r)  =  1  h (r^r1)  $n  dV1 


where 
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h/r  ri\  _  2k2GQ(rtri)U0(ri) 

[L’L  ’ 

G»(r.r‘)  ,G°(|r-r‘|)  .  £g«13jSil> 

4tt  |  J2.-J2 1 1 


The  approximation 


dn  =  (1+  6n)^-l  =  26n  +  6 rr^  ~  26n 


has  been  used. 

I' 

S  The  following  development  (for  a  special  geometry)  is  from  Yeh  and  Lui  (1972) 


region  of 
f luctuat ions 


rece iver 


For  a  spherical  wave 


U0(r)  =  U0(r) 


A0exp(ixr) 

r 


where  r  =  j  r.| 


od 


Mr.r1)  =^-^1  expfik (R+r1  -  r )]. 


From  this,  obtain  the  following  general  expression  for 


Bx  3  <x(ri)x(£*)>  =  k*  . 

/  /B„(r‘  -jj)  *ri  -.r.>] 

V;  vx  n -1  -ul  r  j  Rx 

where  Bn(£{  -  rj)  is  the  covariance  function  for  the  index  of  refraction  fluctua¬ 
tions  which  is  assumed  to  be  a  function  of  coordinate  differences  only.  Note  that 
and  prefer  to  coordinates  at  the  observation  points.  A  similar  expression  for 
Bs  can  be  derived.  This  general  expression  can  be  greatly  simplified  for  specific 
geometries. 


cos[k(R2+  r£-  r*)] 


ri  R., 


dV}  dv£ 
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and  for  the  “forward  scatter  approximation"  wherein  the  scale  size  of  the  fluctua 
tions  is  much  greater  than  A.  Thus  consider  the  slab  geometry  shown  in  figure  2. 


Figure  2. 


transmitter 


a  z 


/.  layer  of  turbulent 

/ ' ' 

•/  fluctuations 


rece iver 


From  Yeh  and  Liu  ( D 97 2 )  B  and  B  can  be  found  from  the  tollow- 

X  s 

i  ng  equat ions : 


BX  (d,0,L)  =  ;  Bs(d;L) 


Ii+Iz 


where 


a+b  b 

1 1  =  /  /  / 
Z=a  z=-b  y=-‘ 


/  Mx,y,z) 


sin  ry'*  ( x*FdZ/L )  ]  dx  dy  d z  dZ 
" 2z( 1-2Z/L ) 


a+b  b 


1-a  z=-b  y-- _  x--t 


Bn'(x,y,z) 

D  +  zz 
T 


sin  [111 

r\  .  Z  J 


dx  dy  dz  dZ 


where  D  =  4Z(L-Z) 
L 


This  expression  for  Bx  (and  Bs)  can  be  used  for  numerical  computations  for  given 
Bn  ,  or  for  simple  forms  of  Bn  and  further  simplifying  assumptions  some  analytical 
expressions  for  the  integrals  may  be  obtained. 

A  similar  geometry  which  allows  for  an  angle  of  incidence  Y  of  the  (unperturb¬ 
ed)  propagation  direction  w.r.t  the  slab  normal  is  used  in  Taylor  (1975)  wherein 
relationships  betv/een  the  spectral  density  functions  for  x  and  n  are  derived  under 
various  simplifying  assumptions  and  for  incident  plane  waves. 
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Thus 


d 

*x(jc)  =  2u  k2secY  /  f(z)  sin2rLj-z  (K^sec2 Y  +  <^)]  . 

^n(Kx.Ky.<xt3nY)dz 

where  f(z)  is  a  function  which  takes  account  of  a  change  of  turbulence  "levels" 
along  the  main  propagation  path.  If  the  layer  slant  thickness,  d,  is  such  that  d  « 
kl2  and  d  «  L  where  1  is  a  mean  scale  size  of  the  fluctuations,  then 


*x[K)  ~  sin2[L  (^xSec^Y  +  *s(£) 

—  2k 

where  ^(k)  is  the  2-dimensional  phase  fluctuation  spectrum  of  the  wave  at  the 
bottom  of  the  slab.  The  same  expression  for  intensity  fluctuation  spectra  is  deriv¬ 
ed  in  the  next  section  indicating  that  under  the  assumptions  used  in  the  derivations 
log-ampl itude  spectra  and  intensity  spectra  are  proportional.  The  function  in 
braces  is  frequently  referred  to  as  the  filter  function.  For  the  thin  layer,  this 
filter  function  is  called  the  Fresnel  filter. 


Section  IV  -  The  thin  phase  screen  approach 

Starting  with  the  parabolic  equation  (2-6),  it  is  possible  to  make  appropriate 
simplifications  (for  various  geometries)  and  derive  expressions  for  the  spectrum  of 
intensity  fluctuations.  In  particular,  the  "thin  phase-screen"  or  simply  "thin 
screen"  approximation  can  be  derived.  The  net  result  is  that  the  spectrum  of  elec¬ 
tron  density  fluctuations  (or  equivalently,  the  spectrum  of  refractive  index  fluctua¬ 
tions)  is  multiplied  by  a  filter  function.  The  derivation  of  the  sin2  filter  was 
first  done  by  Bowhill  (1961)  and  later  by  Salpeter  (1967)  with  results  from  Mercier 
(1962).  The  following  is  an  outline  of  the  "canonical"  thin  screen  done  essentially 
as  in  Ishimaru  (1978). 

Start  with  the  parabolic  equation  (2-6)  for  W(£)  where  the  field  is  written  as 
U(r)  =  W(r)eikz  for  forward  propagation  primarily  in  the  z-direction.  The  intensity 
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fluctuations  are  derived  from  the  fourth  order  moments  of  W(_r) : 


1*4  =  <W(z,Pi)W(z,P2)W*(z,p{)W*(z,P^)> 


With  an  appropriate  change  of  coordinates  from  the  Pi,  p^,  pj ,  p/  system  to  R_,  £, 
n,  n  system,  the  parabolic  equation  for  f4  can  be  put  into  the  form: 


&  -  Q]  '■'>  ■  0 


where 


Q(V^)  =  CH(si)+H(s.J-1/2H(si+sJ-1/2H(si-^)] 

H{fi.)  =  Bn(o)  -  Bn(£). 

£  -  P/£|+P^) 

£  =  £1+£2‘£i“£a 

Sj  =  l^lpj-p  V-pj  ) 

£2  =  l/ZC^-pz-pJ+p^). 

Note  that  the  implicit  assumption  has  been  made  that  only  refractive  index  fluctua¬ 
tions  in  the  plane  perpendicular  to  the  propagation  direction  are  important  in  pro¬ 
ducing  scintillation;  thus  Bn  is  a  function  of  p  =  (x,y)  only. 

The  intensity  fluctuations  near  a  point  z  along  the  propagation  path  j_s 
T  4(z,sJ  ,0). 

The  derivation  of  the  expression  for  F4  at  a  distance  L  from  the  front  of  a 
thin  screen  is  done  in  two  steps  (plane  waves  are  assumed  -  a  correction  for  spheri 
cal  waves  is  given  later).  First,  the  parabolic  equation  is  integrated  across  the 
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thickness  (az)  of  the  screen.  To  first  order,  the  result  is  that  at  the  "bottom"  of 
the  screen  (z= 0) 


r4(0,si,s*)  =  exp(-Q(s.i,s.*),Az) 


(An  alternative  derivation  for  T4  is  to  use  the  Rytov  approximation  to  obtain  an 
expression  for  1’4  at  the  "scintillation  exit  plane  See,  for  example,  Tatarski 
(1961).) 

The  expression  for  T4  at  the  receiver  (z=L)  can  be  obtained  in  various  ways. 
The  result  is: 


r4(L,s1,s_/)  =  (  ><  )z  f  dsjdsj  expdJL(si-sj)  •  (s/-sj)  +  AQ(sj  ,s^> 
L  “oo  L 


Alternatively,  T 4 ( L , s_i ,  _$2 )  can  be  written  as 


oo 

r4(L,si,sJ  =  /  S(L ,£,_s2)e  —  *-A  d< 

“oo 


where 


S(L,k,sJ 


-iK-si 
e - 1 


Note  that  S(L,^,0)  =  ^(L.s^.O)  is  the  spectral  density  function  for  intensity 
scintillations  and  is  essentially  the  quantity  which  is  measured. 

By  making  further  approximations,  it  is  possible  to  derive  various  "filter  func¬ 
tions"  which  associate  the  scintillation  spectrum  at  the  observation  point  (z=L)  to 
that  at  the  phase  screen.  For  example,  if  the  first  order  approximation  eQ  ~  1+Q  as 
used  by  Salpeter  (1967)  is  made,  the  sin*  filter  function  is  found.  Thus: 
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$(L,k,<i>)~  1  c  J  [1  +  AQ(si,-J5L)]e“i- -i  diLi  = 
-  TTi)  -  r 


-iK-s1  .  i 


1  ,  J  [l  -  H(-*L)je  -“i  dIi 
CZ»)  “1“  k 


k2  2AJ  CH(s.[)-l/2H(s.l-^1)-l/2H(s.1+^L)]e“1-‘-5-l  dsj 


The  first  term  gives  rise  to  a  “delta"  function  which  is  unimportant  for  the  study 
of  scintillation  spectra.  The  second  term  can  he  written  as: 

£  - - - 


AH(<)sin^(l^|2L) 


T  ~ 


"IT 


where  [k|2  *  x*  +  y2 

and  H^)  is  the  Fourier  transform  of  H  and  is  essentially  the  (2-dimensional) 
spectrum  of  the  refractive  index  fluctuations.  Thus: 


Si(L,<)  =  S(L,<,o)  ■  C(k)SN(<)sin2(|Jf^L) 


where  C  is  a  constant  which  depends  on  the  free  space  wave  number  k  and  some  geome¬ 
tric  quantities. 

The  above  derivation  assumes  a  plane  wave  incident  on  the  fluctuation  slab  and, 
in  addition,  that  the  propagation  path  is  normal  to  the  slab.  Two  corrections  must 
be  made  for  spherical  waves  and/or  non-normality. 
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(a)  spherical  wave  correction  -  Ishimaru  (1978) 

_ >  +l2  where  Li  is  the  distance  from  source  to  slab  and  L2  is  the 

tT  — 

distance  from  slab  to  receiver. 

(b)  the  angle  between  the  propagation  path  and  slab  normal  is  y.  The  propagation 
path  is  in  the  x-z  plane  -  Taylor  (1975).  kx  — >  k*  secy 

L  — >  slant  distance. 


Section  V  -  A  choice  of  approaches 

As  shown  above,  the  Rytov  and  thin  screen  approximations  give  basically  the 
same  results  for  a  thin  slab  of  fluctuations.  There  is  good  evidence  that  the  thin 
slab  approach  is  applicable  for  radio  waves  traversing  the  ionosphere  wherein  the 
slab  is  located  near  the  electron  density  peak.  In  any  case,  a  choice  has  to  be 
made  and  at  least  for  a  first  approximation  t o  the  truth,  it  is  used  here.  Specifi¬ 
cally,  the  results  derived  in  Taylor  (1975)  for  a  uniform  slab  are  used: 

let  *c  =KxSecY  + 

then 


*(<)  =  irk2  dsec y  {1  +  JS —  Csin(L^d  .  .  sin(L.  *  k 2C  ) }>  <5^ (k x ,  <y xt a n y) 

x  <*c 

For  equitorical  sci nti 1 lation,  $ n  is  highly  anisotropic.  Taylor  (1975)  derives  some 
useful  expressions  under  the  assumption  that  the  correlation  function  for  anisotro¬ 
pic  fluctuations  is  obtained  directly  from  that  for  isotropic  fluctuations  by  a 
scale  multiplication.  The  final  results  are: 
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=  a*n(a2[Kx(cos<(>cose  -  sinetany)  +  <ySin<t)Cos8]  2  +  (-Kxsin4>+KyCos<f>f 
+  [Kx(cos4>sin6  +  cosQtan  ^  +  KySin<J>sirfl  ]2 }  V2  ] 
where: 

^n(<)  is  the  spectral  density  function  for  the  "associated"  isotropic  fluctuation. 

a  is  the  "elongation  factor"  for  the  anisotropic  fluctuation  correl ation  func¬ 
tion.  (Typically,  imagine  correlation  ellipsoids  of  revolution.  Then  a  is 
the  ratio  of  the  major  and  minor  axes.) 

4)  is  the  angle  the  major  axis  plane  makes  with  the  observation  plane. 

0  is  the  angle  the  major  axis  makes  with  the  x,  y  plane. 

The  last  step  is  the  relationship  between (k) ,  the  spatial  spectral  density 
function  and  S^ta),  the  temporal  (observed)  spectral  density  function.  Taylor's 
(1938)  hypothesis  of  "frozen"  turbulence  is  used  to  make  the  transition  (see 
Wyngaard  and  Clifford  (1977)  for  a  discussion  of  the  validity  of  this  hypothesis). 
The  approximation  amounts  to  performing  the  steps: 

(1)  F(KX)  =  fZ  4>(K)dKv,dKz  (for  the  general  3D  case). 

(2)  =  1  F(w)  where  the  "irregularities"  drift  in  the  x-direction  with  speed  v. 
This  approximation  assumes  that  the  x-direction  can  be  freely  chosen.  However, 

for  the  geometry  considered  here,  the  coordinate  system  has  been  fixed  by  the  verti¬ 
cal  plane  through  the  slant  line  joining  transmitter  and  receiver  and  the  local  ver¬ 
tical  . 

Thus  the  following  transformations  must  be  made: 
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(1) <J>  Ux*y)  — >  4>  UxVy)  by  the  substitution 

X  X 

l  1  . 

kx  =  k  xC0S  a  -  KySina 

<y  =  <  xsin  a  +  <ycosoi 

where  the  drift  direction  is  along  the  x1  axis  which  makes  an  angle  a  with  the 
x-axis. 

(2)  F(>ci)  =  /  4>  (icx,Ky)dic^ 

(3)  S(u)  =  i  F(tt) 

V  V 


Section  VI  -  Examples 

The  model  presented  above  contains  many  parameters  which  are  listed  below. 

(1)  general  parameters  which  define  the  isotropic  refractive  index  spectral  den¬ 
sity  function.  For  example,  10  and  L0  for  the  Kolmogorov  spectrum. 

(2)  d,  the  slab  slant  thickness. 

(3)  L,  the  slant  range  to  the  top  of  the  slab  from  the  receiver. 

(4)  the  velocity  of  the  slab  relative  to  the  slant  mean  propagation  path.  This 
consists  of  two  parameters,  v  and  a  . 

(5)  the  parameters  a,e  and  which  describe  the  anisotropy  of  the  index  of 
refraction  fluctuations  relative  to  the  coordinate  system  defined  by  the 
slant  path. 

The  coordinate  system  defined  by  the  slant  path  changes  in  time  due  to  motion 
of  the  transmitter  and  receiver.  In  addition,  of  course,  the  fluctuation  slab  is 
assumed  to  move  relative  to  the  slant  path.  (The  condition  v=o  can  presumably  oc¬ 
cur,  in  which  case,  the  scintillation  spectrum  S  is  not  defined.)  During  a  given 
observation  of  S,  it  is  assumed  that  all  parameters  remain  constant.  A  computer  sub¬ 
routine  has  been  written  (see  Appendix  A)  which  provides  the  parameters  v  and  a 
given  the  satellite  (orbital  parameters  X,Y,Z,  X,Y,2),  the  aircraft  position  and 
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velocity  vector  for  the  mean  slab  motion  (relative  to  a  fixed  earth).  Various 
spectra  have  been  computed  using  the  above  models  with  various  values  for  the  para¬ 
meters.  In  particular,  these  calculations  have  been  done  for  Gaussian,  power  law 
and  Kolmogorov  (or  von  Karman)  spectra. 

Figures  3a,  b,  c  show  the  function  o>x{<x,<y)  for 

a)  Gaussian 

b)  vn  power  law 

c)  in  von  Karman 

The  following  parameter  values  for  4>x  were  used: 


transmitter  frequency 

100MHz 

slant  angle  (y) 

OO 

distance  from  layer 
to  receiver  (L) 

300  km 

layer  thickness  (D) 

10km 

correlation  "ellipsoid" 

values : 

(a) 

2. 

(e) 

2b  deg. 

(♦) 

45  deg. 
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It  should  be  noted  that  all  of  the  physics  of  this  particular  model  of  ampli¬ 
tude  scintillation  is  contained  in  these  functions,  (kx,kv).  The  amplitude 
spectra  are  obtained  from  these  functions  by  first  rotating  the  <-coordinate  frame 
over  which  they  are  defined  by  angle  a  and  then  integrating  from  -00  to  00  the  result¬ 
ing  function  in  the  <y  direction  for  all  values  of  kx.  Because  of  symmetry,  values 
of  Kxl  2  0.  only,  need  be  considered.  The  final  "transformation"  from  *  (Kx»*y)  to 
S(w)  is  merely  one  of  stretching  or  compression  and  normalization  such  that  the  area 
under  the  function  remains  the  same.  Thus  it  is  important  to  study  the  structure  of 
these  functions  in  both  a  qualitative  and  quantitative  manner.  The  transformations 
mentioned  above  can  be  done  mentally  in  a  qualitative  fashion  using  the  illustrated 
functions  as  examples  of  In  this  report,  the  structure  of  $xis  not  investi¬ 

gated. 


In  Figure  4,  which  is  a  plot  of  the  function  S(f)  =  S(^),  some  examples  are 
shown  of  the  above  transformations  of  a  single  for  an  actual  aircraft/satellite 
geometry  (as  occurred  on  October  20,  1976).  This  figure  is  an  illustration  of  the 
output  of  a  set  of  computer  software  designed  for  comparing  theoretical  and  experi¬ 
mental  spectra.  Again,  it  is  not  the  intent  of  this  report  to  delve  into  the  struc¬ 
ture  of  the  spectra.  Appendices  A  and  B  provide  some  information  on  the  algorithms 
used  in  producing  Figure  4. 

The  following  fixed  parameters  were  used: 


satel 1 i te 
transmitter  freq. 
date/time 
A/C  height 
A/C  location 
A/C  speed 
layer  model 
layer  speed 
layer  di recti  on 
layer  "shape" 
layer  thickness 


LES9 
249  Wz 

76/10/20/1/8/0 
30000  ft. 

lat  -lio,  long  -77° 

470  knots 
power  law  (p=3) 

70  m/s 

900  (geographic  east) 

9  =  0°,  <J>  =  0°,  a  =  10 

10  km 
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Section  VII 


In  this  section  a  brief  overview  of  "spectral  analysis"  is  presented  along  with 
a  description  of  the  actual  spectral  analysis  technique  used  to  compare  theory  with 
experiment. 

The  theoretical  quantity  of  interest  is  the  "measurable"  scintillation  spectral 
density  function,  S(w),  related  as  shown  previously  to  the  spectral  density  function 
of  electron  density  fluctuations.  The  function  S(co)  is  presumably  well  defined  and 
measurable  (or  more  exactly  determinable)  over  a  portion  of  the  frequency  domain  - » 
to  +co.  Because  the  underlying  stochastic  process  is  real,  S(oj)  is  even  and  need 
be  determined  only  over  o  to  ».  Moreover,  S(a$,  for  all  practical  purposes,  is  negli¬ 
gible  beyond  some  to  which  can  be  determined  from  the  model  parameters  as  listed  in 
the  previous  section.  S(w)  is,  of  course,  a  probabilistic  concept,  summarizing  some 
(and  perhaps  all)  of  the  second  moment  properties  of  an  underlying  stationary  sto¬ 
chastic  process.  The  task  at  hand  is  to  estimate  S(^)  or  perhaps  a  portion  of  S(oj) 
given  a  finite  (and  relatively  small)  observation  (or  "realization")  of  the  process. 

Thus  we  enter  the  broad  expanse  of  another  discipline  called  spectrum  (or  spec¬ 
tral)  analysis  or  more  exactly  (but  verbose),  the  theory  and  practice  of  estimating 
a  spectral  density  function  of  a  (in  general)  "homogeneous  process"  given  one  or 
more  finite  realizations  of  this  process.  The  term  "homogeneous  process"  is  put  in 
quotes  because  quite  often  it  is  not  clear  that  the  underlying  model  is  a  homogene¬ 
ous  process  (e.g.,  to  what  extent  is  human  speech  a  homogeneous  process),  and  modern 
spectral  analysis  is  not  limited  to  homogeneous  (or  stationary)  processes  although 
it  often  is  not  clear  how  the  spectrum  should  be  defined  for  "non-homogeneous"  pro¬ 
cesses.  None-the-less,  it  is  quite  clear  in  this  study  of  scintillation  spectral 
analysis,  the  process  behind  the  S(oj)  is  (assumed  to  be)  stationary.  (Note  that  the 
term  homogeneous  usually  is  reserved  for  stochastic  processes  defined  in  space  while 
stationary  is  reserved  for  "random  signals",  i.e.,  "random  functions  of  time".) 

The  underlying  process  (random  signal  or  whatever  one  wants  to  call  it)  is  an 
observed  AVC  voltage  which  can,  by  suitable  transformations,  be  made  proportional  to 
the  scintillation  amplitude  A  as  defined  in  the  glossary.  Note  that  synchronous 
detection  is  not  used  here  so  that  the  X  and  Y  components  of  the  signal  cannot  be 
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determined  separately.  In  addition  to  being  of  finite  duration,  the  process  is  digi 
tized  either  directly  during  the  recording  process  (after  Jan.  1979)  or  prior  to 
analysis.  (Prior  to  1979,  the  signal  voltage  was  recorded  directly  on  an  FM  magnetic 
tape  recorder  and  later  digitized  prior  to  analysis).  Thus  the  underlying  process 
can  better  be  thought  of  as  a  "time  series"  and  the  usual  techniques  of  time  series 
spectrum  analysis  can  be  applied. 

Before  delving  into  the  subject  of  time  series  spectrum  analysis,  it  is  worth¬ 
while  to  bridge  the  gap  between  the  "continuous  domain"  and  the  "discrete  domain" 
brought  about  by  sampling  and  to  then  forget  entirely  about  the  continuous  domain 
(it  remains  tied  to  the  discrete  domain  through  some  frequency  conversion  constants 
and  assumptions  on  the  "band  width"  of  the  original  signal). 

First  of  all  it  is  assumed  that  the  sample  interval  At  (assumed  constant)  is 
always  such  that  the  Nyquist  criteria  is  met,  i.e.,  At  £  tt/wc  where  wc  is  such  that 
S(oj)  =  o  for  j^i  >  In  practice,  this  condition  is  better  stated  as  "the 
sample  interval  is  small  enough  such  that  spectral  aliasing  is  negligible."  Unfor¬ 
tunately  this  condition  may  sometimes  be  inadvertantly  violated.  For  example,  un¬ 
known  high  frequency  noise  may  be  present  in  the  observing  system.  Violation  of 
this  condition  can  frequently  be  detected  by  estimating  the  spectrum  at  various 
sample  intervals. 


Once  in  the  "discrete  domain"  it  is  convenient  to  use  a  unit  sample  interval 
(basically  a  time  series  is  a  random  function  over  the  integers  (usually  represented 
in  group  theory,  at  least, by  Z)).  The  spectral  density  functions  which  might  be 
defined  for  such  time  series  is  defined  over  -  tt  <  o>  <  ^  (more  exactly  S (a> )  is 
periodic  with  period  2tr).  The  estimated  spectral  density  functions  which  are  defin¬ 
ed  over[-Tr,Tr  ]  can  be  associated  with  the  original  (unsampled)  process  by  using 
the  conversion  factor^.  ,  i.e.,  S  (  oj)->S  (^t") 

In  particular,  given  the  sample  interval  At,  the  spectrum  can  only  be  estimated  up 
to  a  (radial)  frequency  given  by 

jY  (or  up  to  ^  Hertz)*  It  is  intuitively  obvious  perhaps  that  the  sample  size,  N, 

(the  original  observation  now  consists  of  N  signal  values  x0,  xx .  xn_x  )  limits 

the  low  end  of  the  spectrum  or,  more  exactly,  essentially  limits  the  accuracy  of 
estimating  the  spectrum  nearoj  =  0.  [This  is  usually  not  a  problem  since  one  usually 


is  not  interested  nor  should  be  interested  in  S(w)  at  u>=  0  .]  Roughly  speaking,  the 
lower  limiting  (radial)  frequency  for  a  reliable  estimate  of  S(u>)  is  w  where 

M  <  N  is  a  "truncation"  number  used  to  reduce  the  variance  in  the  estimate  of  S. 
(Because  of  the  previously  mentioned  evenness  of  S(w)  it  is  only  necessary  to  dis¬ 
cuss  estimating  S  for  positive  frequencies.) 

Historically,  spectral  analysis  of  time  series  really  got  its  main  "modern" 
impetus  with  the  publication  in  1959  of  "The  Measurement  of  Power  Spectra",  Blackman 
and  Tukey  (1959).  During  the  decade  of  the  60's  the  theory  and  technique  of  spec¬ 
tral  analysis  became  the  subject  of  many  articles  and  books.  During  this  period, 
the  subject  arrived  at  its  "classical  formulation"  which  is  still  the  one  most  often 
used.  (It  is  quite  adequate  for  many  purposes.)  A  good  summary  of  this  formulation 
can  be  found  in  Koopmans  (1974). 

During  the  last  decade,  more  "powerful"  techniques  of  spectral  analysis  have 
become  more  common  place.  These  methods  are  usually  referred  to  as  the  MEM  or  MLM 
methods  (for  maximum  entropy  method  and  maximum  likelihood  method,  respectively). 
These  (currently  modern)  methods  as  well  as  the  "classical"  methods  of  spectral  anal¬ 
ysis  are  summarized  in  a  collection  of  reprints  by  Childers  (1978).  The  Proceedings 
of  the  RADC  Spectrum  Estimation  Workshop  (1978)  also  contains  several  up-to-date 
"working"  methods  for  spectrum  analysis  and  provides  an  overview  of  many  diverse 
areas  where  spectral  analysis  is  used.  A  fairly  recent  article  which  gives  the 
"modern  flavor"  of  spectral  analysis  can  be  found  in  Thomson  (1977). 

Before  giving  the  particular  algorithm  which  is  used  here  for  estimating  S(w) 
from  the  observed  time  series,  it  is  worthwhile  to  first  quickly  review  some  of  the 
primary  ideas  which  led  to  the  development  of  the  algorithm  which,  by  the  way,  is 
approximately  equivalent  to  most  other  "classical"  methods. 

First  of  all,  the  whole  idea  of  associating  a  "spectrum"  with  a  time  series  is 
an  attempt  to  "Fourier  analyze"  (or  harmonic  analyze)  observed  data  -  to  break  it 
down  into  its  fundamental  components,  so  to  speak.  Mathematically,  the  problem  is 
one  of  generalizing  the  concept  of  Fourier  series  or  the  Fourier  integral  for  finite 
"energy  functions".  This  generalization  can  be  approached  from  two  directions, 
i.e.,  through  probability  theory  (stationary  stochastic  processes)  or  finite  "power" 
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but  "deterministic"  functions.  (From  a  practical  standpoint,  there  is  little  obvi¬ 
ous  difference  between  a  (segment  of)  a  deterministic  waveform  and  a  (segment  of)  a 
realization  of  a  stochastic  waveform.  The  main  difference  is  that  with  the  latter 
one  can  associate  various  probabilities  and  "expected  values",  while  with  the  former 
this  cannot  be  done,  at  least  not  with  the  same  interpretations.)  The  approach  to 
spectral  analysis  via  stochastic  process  theory  is  less  intuitive  than  that  through 
finite  power  function  theory  so  we  choose  the  latter  to  start  with  but  soon  switch 
over  to  results  from  stochastic  process  theory.  This  approach  also  has  the  advan¬ 
tage  that  the  rationale  behind  the  spectrum  estimation  algorithm  is  somewhat  more 
obvious. 

To  start  with,  consider  the  class  of  real  or  complex  valued  functions  on  Z  (the 
integers)  such  that 


(7.1)  <p  (n)  = 


nm 


1 


N-*-00  2N+1 


N 
E 

n'=-N 


f(n')f(n*+n) 


exist  for  all  n  e  Z. 

(This  class  is  the  discrete  analog  of  "finite  power"  signals  f(t)  such  that 


<D(t)  = 


1  im 

J-+-0O 


1 

7T 


f  (t ' )  f  ( t '  +t )  dt ' 


exists  and  is  non-trivial  for  all  t  e  R.  Since  we  deal  only  with  "digitized"  wave 
forms,  it  is  simpler  to  consider  functions  defined  over  Z  to  begin  with  but  to  keep 
in  mind  the  sample  interval  1  tt  and  the  specter  of  possible  aliasing.)  4(n)  is 
called  the  "autocorrelation  function  of  f".  Functions  of  this  class  were  first 
seriously  considered  by  Norbert  Wiener.  They  are  not  amenable  to  harmonic  analysis 
in  the  usual  sense  although  it  is  possible  to  develop  a  "generalized"  harmonic  analy¬ 
sis  for  such  functions  which  is  very  similar  to  that  which  has  been  developed  for 
stochastic  processes.  The  function  <j(n)  is  positive  definite  and  can  always  be  rep¬ 
resented  in  the  form 
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4>(n)  =  eina)dF(to)  and  if  F  is  absolutely  continuous, 
<|>(n)  =  fn  s(u>)  e’nw  dw.  Furthermore,  assuming  the  inverse 

-IT 

Fourier  transform  exists. 


(7.2)  s(w)  =  ^  ^  4>(n)-in<^ . 


The  function  s(ui)  is  frequently  called  the  spectral  density  function  of  f;  it 
is  the  functional  analog  of  S(w),  the  spectral  density  function  associated  with  a 
stationary  stochastic  sequence  (or  time  series). 

Non-trivial  examples  of  f  with  a  given  associated  s  can  be  constructed  using 
"Wiener  numbers",  see  Papoulis  (1962)  as  input  to  a  digital  filter  with  amplitude 
response  |H(oj)  j  such  that  | H (oj )  |  -  s(u>). 

Consider  next  the  problem  of  determining  (or  estimating  s(oi))  given  f(n)  on 
[0,N]  or  [-N,N].  More  explicitly,  consider  a  sequence  of  (in  this  case,  real 
valued)  functions  f|y(n)  on  Z  which  are  identically  zero  for  |  N|  >N  and  such  that 
fN(n)  =  f(n)  f°r  (nliN*  For  each  such  fj<  consider  the  quantity 


(7.3)  SNU  71PT  !  fN(n)e-inu>!2. 


Using  results  from  Fourier  analysis,  it  is  easy  to  show  that 


(7.4) 


SNk>  ) 


2N 

Z  ♦N(n)e"inu> 
-2N 


where 


40 


X 


4>N(n)  =  TTJ+r  ^N*fN(n)  (where  *  =  correlation) 


i.e.,  in  general,  f*g(n)  =  £  Ffn'Jgfn+n1)  for 

n'eZ 

complex  valued  f,g  and  ~  means  conjugation. 


Comparing  equations  (7.4)  and  (7.2),  it  seems  plausible  that  s(w)  can  be 
"estimated''by 


11  * 

(7.5)  |  E  fN(n)e“'irw 1 2  for  some  suitably  large  N. 

n=-N  ! 


We  now  turn  immediately  to  the  completely  analogous  situation  where  instead  of  a 
(deterministic)  function  f  ,  we  have  a  time  series  (realization  of  a  stationary  sto¬ 
chastic  sequence),  specifically  N  values  of  the  time  series  and  define  the  analog 
of  equation  (7.5). 


(7.6)  InH  =  N  |  n£1XN(n)e-ina)| 


(Here  the  time  series  x  is,  for  historical  reasons,  indexed  from  1  to  N  rather  than 
from  -N  to  N.) 


h  £  RN(n)e-->"“ 

n— N+l 


1  1  N_1 

where  Rn( n )  =  Rm(-n)  =  XN*XN(n)  =  w  z  XN(n‘ )xn( n ’ +n ) 

M  n’  =  l 


N-n 

Z  x(n')x(n'+n) 
N '  =1 
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where  xn(h)  =  x(n),  n=l,..., 


N;  xn(ii)=0  elsewhere. 


In  is  called  the  "periodogram"  and  Rn  the  correloqram  associated  with  the  time 
series.  They  are  the  basis  for  most  of  “classical"  spectral  analysis. 

It  can  be  shown,  see  for  example  Hannan  (1967),  that  In  is  not  a  "consistent" 
estimator  of  S(gj)  associated  with  the  time  series.  That  is,  the  variance  of  In  does 
not  approach  zero  as  the  length,  N,  of  the  time  series  goes  to«>  (although  the  expect¬ 
ed  value  of  In(u)  approaches  S(w)).  In  fact,  at  least  for  a  Gaussian  time  series 

var[INfo  )]  - >  a4  S2  ( co) 

where a2 is  the  variance  of  x.  (Unless  noted  otherwise,  it  is  assumed  that  x  has 
mean  zero.) 

Faced  with  this  "inconsistency"  dilemma,  the  most  obvious  solution  is  to  some¬ 
how  obtain  an  "average"  periodogram  in  the  hope  that  the  "fluctuations"  in  the  esti¬ 
mator  (as  measured  by  the  variance  of  the  estimator)  will  approach  zero  for  large  N. 
In  fact,  the  art  of  spectral  analysis  (specifically  the  art  of  spectral  density  func¬ 
tion  estimation)  was  initiated  by  Bartlett  (1950)  when  he  published  results  of  a 
careful  analysis  of  this  averaging  process.  Briefly,  the  averaging  is  done  by  break¬ 
ing  the  time  series  of  length  N  into  M  segments  of  equal  length  N1,  compute  the 
periodogram  for  each  segment,  and  then  average  them.  Without  going  into  the  details, 
it  can  be  shown  that  this  process  is  equivalent  to  estimating  S(cj)  by 

R  1  NM 

sN(w)  =  Z  WN(n)RN(n)e"in w 

n  n=-N'+l 

where  wn ( n )  =  1-  j  n  | /N '  -N‘  _<  n  ^  N' 

=  0  elsewhere 
and  N  =  MN' 

B  B 

The  variance  of  Sn,  (var[SN(u> )] ,  is  given  by 


var[s[j(  id)]  =  var[If^(a))]/M. 

Thus  assuming  that  M  becomes  large  with  N,  is  a  consistent  estimator.  The  next 
logical  step  is  to  introduce  a  general  estimator  of  the  form 


(7.7)  Sn(u>)  =  jz  I  w(n)R/v(n)e"1n“ 

n— mt+1 


4-  £  WN(n)RN(n)e-'' 

a  neZ 


where  wn  is  defined  to  be  zero  for  | n  | >N .  The  integer  mt  is  referred  to  as  the 
"truncation  point"  and  the  function  wn  is  called  the  "lag  window". 

A  further  generalization  is  developed  in  Grenander  and  Rosenblatt  (1957)  where 
the  estimator 


1 

snU)  l  t>N(n,n')x(n)x(n')  is  analyzed. 

n,n‘ 

The  bN(n,n‘)  =  bN(n',n)  is  a  set  (indexed  by  N)  of  real  valued  functions  on  Z  x 
Z.  Furthermore,  if  bN(n,n')  is  chosen  such  that 


bN(n.n')  =  V^yje1  (n“n '  )ydy 


where  Vn  is  a  real  or  complex  valued  function  with  a  Fourier  transform,  then  one 
obtains  essentially  an  estimate  of  the  form  given  by  equation  (7.6).  Estimates  of 
this  type  are  called  spectrogram  estimates.  If  vfj  (the  Fourier  transform  of  Vn)  is 
identically  zero  for  |n(>my-l,  then  Sn(uj)  is  called  a  truncated  spectrograph  esti¬ 
mate. 
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One  advantage  of  spectrograph  estimates  (as  a  class)  is  that  their  statistical 
properties  are  quite  well  understood.  There  are  several  other  "classical"  methods 
of  spectral  analysis  which  are  not  strictly  of  the  spectrograph  type.  For  example, 
a  common  method  is  to  divide  the  time  series  of  length  N  into  M  segments  which  gener 
ally  overlap.  A  suitable  real  weight  function  (or  "fader"),  a(n),  is  chosen  and  the 
quanti ty 

2  m- 1  .2M 

zi (j)  =  l  a(n)x] (n)e1  2m 
n=0 


j=0,...m  is  found  for  each  segment  1  where  2m  is  the  length  of  each  segment. 
The  quantities  |z](j)  |  2  are  then  averaged  to  provide  estimates  of  Sn(w)  for 


In  this  method  the  data  itself  is  "windowed"  rather  than  the  correlogram. 

It  should  be  pointed  out  here  that  all  of  the  "classical"  methods  are  "ad  hoc" 
in  the  sense  that  an  estimator  (or  estimation  method)  is  proposed  (e.g.,  the  general 
spectrograph  estimator)  and  then  some  of  its  statistical  and  other  properties  are 
investigated.  All  of  them  include  a  "window  function"  of  some  sort  and  much  has 
been  written  about  the  properties  of  various  windows.  Also  these  methods  share  the 
common  practice  of  appending  zeros  to  the  data  (or  the  correlogram),  a  practice 
which  can  be  paraphrased  as  observer  introduced  "bias"  to  the  data.  Arbitrary  win¬ 
dowing  and  adding  zeros  are  both  practices  which  is  tantamount  to  the  user  inserting 
"information"  in  the  estimating  procedure  -  information  which  has  no  basis  in  fact. 
The  introduction  to  the  modern  era  of  spectral  analysis  is  aptly  summarized  by  Abies 
(1978)  wherein  any  data  reduction  method  (including  spectral  analysis)  should  be 
"consistent  with  all  relevant  data  and  maximally  non-committal  with  regard  to  una¬ 
vailable  data"  (roughly  speaking,  a  maximum  entropy  principle).  Unfortunately,  this 
principle  which  is  so  appealing  and  easy  to  state  does  not  provide  for  a  method  of 
spectral  analysis;  one  still  must  devise  a  method  which  "works"  in  some  sense  and 
then  see  if  it  follows  the  maximum  entropy  principle. 


For  the  task  at  hand,  we  do  not  really  want  to  compute  the  spectrum  per  se,  but 
to  do  some  hypothesis  testing,  e.g.,  test  the  hypothesis  that  the  observed  spectrum 
is  due  to  a  power  law  electron  density  spectrum  against  the  alternative  hypothesis 
of  a  Gaussian  electron  density  spectrum,  etc.  Never-the-less ,  in  this  report,  the 
hypothesis  testing  approach  is  not  used  from  "first  principles".  Rather  theoretical 
spectra  and  experimental  spectra  are  derived  and  compared  in  the  obvious  manner. 
Moreover,  for  reasons  of  expediency,  the  classical  (truncated)  spectrograph  estima¬ 
tor  is  used.  Details  of  the  estimator  (estimation  method)  are  given  in  Appendix  C. 

I 

This  section  on  spectral  analysis  is  completed  by  presenting  some  well  known 
properties  of  spectrograph  estimators  (derivations  and  more  details  can  be  found  in 
Koopmans  (1974)). 

(1)  Confidence  region 

If  it  is  possible  to  derive  a  probability  distribution  function  for  a  given 
estimator  Sn(w)  and  to  do  this  for  all  ( 0  <  oi  1  n ) »  then  calculated 
estimates  can  be  placed  within  a  confidence  region  or  more  exactly  a  confidence 
region  can  be  superimposed  on  a  set  of  estimates  Sft(co)  (note  that  a  generic 
estimator  (actually  a  random  function)  is  denoted  by  Sn(J,  while  a  specific 
estimate  is,  here,  denoted  by  Sft(oj)).  It  is  then  possible  to  make  the  pseudo- 
mathematical  statement  that  "with  P%  confidence,  the  true  value  of  S(w)  is  with 
in  the  given  region"  (which  depends  among  other  things  on  S*(u>)  itself). 

Approximately,  it  can  be  shown  (see,  for  example,  Koopmans  (1974)  or  Jenkins 
and  Watts  (1968)  that  the  probability  distribution  for  |^^-is  given  by  P^x^) 
where  P(x|u)  is  the  chi-squared  distribution  for  u-degrees  of  freedom  andu  is 
gi  ven  by 


2N 

U  ’  «T|  M  £ 


(u  is  half  the  above  value  at  the  ends  (o.tt)  of  :he 
i  nterval ) 


where  I  I  h I  I  =  [  / 1  h2 (t )dt] x/2  is  the  norm  of  the  window  function  referred  to 

2  — ^ 

previously.  (It  is  customary  to  define  the  various  window  functions  such  that 
Wf((n)=  h (|f)y  where  h  is  defined  on  R  with  support  [-1,1].) 
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I 


i 


It  is  convenient  to  display  log  S^*(ai)  rather  than  Sn*(w)  since  the  length 
of  the  confidence  interval  for  log  S(w)  is  independent  of  m  .  The  100(l-a)% 
confidence  interval  is 


[logQJ)  +  log  Sn*(oj),  log(u)  +  1  ogSN*(o>)] 


where  a  and  b  are  such  that 


P(ux  |  u  <  a)  =  %-> 

P  (u  x  |  u  b )  =  1  - 

From  this  it  follows  that  log  &  is  the  constant  length  of  the  confidence 

a 

interval  on  a  log  scale. 

In  deriving  the  above  results,  several  assumptions  are  made  but  which  are 
not  strictly  true.  For  example,  it  is  assumed  that  the  bias  of  the  estimator 
(see  next  sub-section)  is  zero;  it  is  also  assumed  that  the  estimates  at  the 
various  w  values  are  uncorrelated  (also  see  below).  Actually  estimates  are 
approximately  uncorrelated  only  at 


“  =wk  =  for  k  =  0,...  mj/2. 

(2)  Bias 

The  bias  of  an  estimator  is  a  measure  of  how  "closely"  the  estimator  can,  in 
the  limit,  estimate.  This  measure  is 


Bn ( a )  =  E  SnGd)  -  S (cjj )  (E  means  expectation) 
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It  is  usually  difficult  to  find  useful  expressions  for  bias  and  they  frequently 
depend  on  derivatives  of  S ( oj)  which  are  unknown  (since  S^)  is  unknown).  For 
example,  for  the  window  function  used  in  this  report 


Bn(w)~  f  S "(m)  ("  means  second  derivative) 

Wj 


A  very  important  consequence  of  bias  (which  usually  can  be  ignored  for  rela¬ 
tively  smooth  spectra)  is  the  inaccuracies  it  produces  when  the  observations 
(the  time  series)  have  non-zero  mean.  In  fact,  if  the  mean  of  the  series  is 
large  compared  with  its  variance,  bias  will  cause  the  estimator  to  behave  like 
(l-cosw)-1  for  all  w  regardless  of  the  true  shape  of  S(w). 

(3)  Variance 


Like  bias,  it  is  difficult  to  obtain  an  "exact"  expression  for  the  variance 
of  an  estimator.  An  approximate  form  (which  follows  from  the  chi-squared  proba¬ 
bility  distribution  for  $n(gj))  is 

var[SN(w)]  ~  —  S2  (w)|  |  h  |  |  *  u>=0,tt 

S  (w)|  | h  |  |  ^  oo^O,  it 


(4)  Correlation 


The  covariance  of  two  random  variables  is  an  important  parameter  since  non¬ 
zero  covariance  at  least  indicates  that  the  random  variables  cannot  be  indepen¬ 
dent.  Covariance  is  frequently  measured  by  the  correlation  coefficient  given 

by 


covCSn^j  ),SN(W2)]/(var[SN(U)i)]  •  var[SN(  w2)])y2 
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where 


cov[x,y]  =  E  [(x-£x)(y-Ey)]. 

Intuitively,  at  least,  covariance  is  a  measure  of  resolution,  i.e.,  if  Sn(uu) 
and  Sn (co 2 )  have  non-zero  correlation,  they  are  somewhat  "contaminated"  by  one 
another.  It  can  he  shown  that  asymptotically  SnUj)  and  SnCjj2  )  are 
uncorrelated  if  | oi i  -  u>2  |  2^:  (independent  of  the  laq  window  used). 

(5)  Bandwidth 


Bandwidth  is  another  measure  (as  is  correlation)  of  the  ability  of  to 
"resolve"  peaks  of  S.  It  is  a  non-statistical  concept  which  is  window  depen¬ 
dent.  There  are  various  measures  of  bandwidth  used  in  spectral  analysis  but 
they  all  are  based  on  the  fact  that  spectrograph  estimator  can,  at  least  approx 
imately,  be  put  in  the  form 

sn(J  -f  wN(y-  )iN(y)dy 

-rr 

where  Wn  is  called  the  spectral  window  (the  Fourier  transform  of  the  lag 
window).  If  one  considers  the  periodogram  In  as  containing  the  maximum  "inform 
ation"  about  S,  then  this  representation  states  that  Sn  is  obtained  by  "smear¬ 
ing"  or  averaging  the  periodogram  over  a  band  of  frequencies  given  roughly  by  a 
width  figure  for  Wn  (and  hence  the  reason  for  its  name  "spectral  window").  A 
frequently  used  width  figure,  g  ,  (see  Parzen  (1961)  is  given  by  1/2  the  base  of 
a  rectangle  with  the  same  area  as  Wn  and  with  height  Wn(o).  For  a  truncated 
estimator  this  is 

f»[SN]  =  ~  ^h(t)dt 


4P 


. . 'mm 


up 


where  h  is  as  defined  above.  Usually  Z1  h(t)dt  is  approximately  1,  thus  it 
can  be  seen  that  bandwidth  is  approximately  equal  to  the  correlation 
"distance",  i.e.,  the  distance  between  uncorrelated  estimates  Sn(wi),  S^(co2). 


It  is  important  to  note  that  resolution  or  correlation  is  related  to  the 
performance  of  Sn(oj)  at  id  =  0.  Strictly  speaking,  Sn(w)  is  statistically  only 
twice  as  "bad"  at  0  and  2tt  than  elsewhere.  (This  factor  of  2  is  due  essen- 

TTTf  ' 

tially  to  the  fact  that,  since  S  is  even,  it  can  be  considered  as  folded  at w = 

0  and  the  estimate  of  S(-w)  is  "folded  into"  the  estimate  of  S(w).)  On  the 
other  hand,  one  can  say  that  the  estimate  S^(e)  where  0  <  e  <  —  is  no  good 

because  of  bandwidth  considerations.  One  can  also  interpret  bandwidth  as  mean¬ 
ing  that  since  one  cannot  "resolve"  frequencies  which  are  closer  together  than 

2ir  2  n  271 

— ,  one  cannot,  in  fact  resolve  any  frequency  less  than  — ,  i.e.,  —  is  the 

lowest  frequency  for  which  one  can  obtain  a  reliable  estimate.  As  a  matter  of 
fact,  this  is  an  optimistic  estimate,  since  usually  the  time  series  under  inves¬ 
tigation  is  not  strictly  "stationary".  For  example,  the  time  series  typically 
has  a  "slowly  varying  mean  value"  (whatever  this  means)  which  must  be  removed 
prior  to  spectral  analysis.  The  removal  process  is  usually  not  completely  suc¬ 
cessful,  i.e.,  there  remains  a  residual  non-zero  mean  which  causes  the  esti¬ 
mates  nearu)=0  to  be  strongly  biased.  Thus  to  be  on  the  safe  side,  the  limit- 
ing  low  frequency  should  be  several  cimes 
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APPENDIX  A 


Calculation  of  the  Parameters  y,  L,  v,  a 


* 


Figure  1 

The  position  and  motion  of  the  satellite  at  point  2  are  defined  by  its 
coordinates  X,Y,Z  in  a  (celestial)  inertial  coordinate  system,  and  its  velocity 
components  X,Y,Z  in  this  system. 

The  position  and  motion  of  the  aircraft  at  point  1  are  defined  by  aircraft 
height  (ha),  geographic  latitude  ($a),  geographic  longitude  (Aa),  ground  speed  (va) 
and  direction  with  respect  to  geographic  north  (aa).  (The  direction  is  assumed  to 
be  the  direction  of  the  ground  path  -  not  necessarily  the  actual  aircraft  heading. 

SC 


The  scintillation  layer  is  "defined11  by  its  height  above  the  earth  hs  and  its 
motion  by  vs,  a$  where  vs  is  its  speed  of  motion  and  as  its  direction  of  motion  with 
respect  to  geographic  north. 

The  point  P  is  the  point  where  the  transmission  path  (lineTT)  intersects  the 
scintillation  layer.  The  motion  of  the  point  P  relative  to  the  scintillation  layer 
must  be  found  as  well  as  the  distance,  L,  between  P  and  1,  and  the  angle, y  ,  between 
the  zenith  at  P  and  the  transmission  path.  The  motion  of  P  is  to  be  obtained  as 
speed,  v,  and  direction, a  ,  relative  to  the  local  x-axis  defined  to  lie  in  the  plane 
defined  by  P  ,2  and  P , zenith. 

It  is  assumed  that  during  an  observation  period,  the  quantities  y  ,  L,  v  and  a 
remain  constant. 

In  the  following  derivations,  an  earth  centered  coordinate  system,  fixed  to  the 
earth,  is  used  (this  is  called  the  computation  system).  The  relationship  between 
this  system  (x,y,z)  and  the  celestial  system  (X,Y,Z)  is  given  by  the  Greenwich 
Sidereal  Time  (here  defined  as  an  angle,  &q)  and  the  earth's  rotational  speed 

1  e  =  It  * 

For  determining  the  aircraft  position  in  the  computation  system,  the  earth  is 
considered  to  be  the  standard  ellipsoid.  All  units  are  assumed  to  be  compatible, 
e.g.,  distances  in  kilometers,  angles  in  radians  and  time  in  seconds. 
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(a)  Conversion  of  satellite  position  and  motion  to  the  computation  system. 


x2  =  Xcos6g  +  Ysinee 

y2  =  -XsinSQ  +  Ycos0g 

z  =  Z 
2 

*2  =  Xcos0g  +  Ysi n@G  +  y2ae 
^2  =  -Xsin0G  +  Ycosog  "  x2  ^ 

K  -  l 

(b)  Conversion  of  aircraft  position  and  notion  to  the  computation  system. 

\  =  (u+ha)cos<j>acosAa 
y%  =  (u  +ha  )cos4)asi  nX  a 
Zj  =  [u  (l-e2)+ha]sin<j>a 

where 

_ ^m _ 

( l-e2s i  n2 0 a ) 1  /  2 

arn  =  semi  major  axis  =  6378.  16  km 
£  =  2f-f2  =  eccentricity  where  f  = 

Xj  =  -(u+ha)[cos  4fcsinAaXa+  s  in4>a<JacosA  a] 

£  =  (u+ha)[cos  (fecos  AaAa-  si  n^si  nAa] 

ii  =  [u(  1-e2  )+ha]cos<{^a 

where 

•  va 

<j,a  ^p-cosoa 

Ka 

.  __  va  sinu, 

Aa~  cosVa 

Ra  =  u  +  ha 

(c)  Computation  of  the  point  P  coordinates  and  motion. 


(2)  "P.P  =  Rp2  ~(^hs)z 

where  ?,  Pi  ,  P2  are  the  vectors  to  points  P,  1,  2  from  the  origin.  The 
quantity  0<  a<  1  is  to  be  computed. 

Substituting  (1)  into  (2),  the  following  expression  for  a  results: 


a  a2  +  b  a  +  c  = 

(R*  +  R?  -2^1  *^)a2  +  2(1*  -R2  )a+R^  -Rp  =  0 


Solve  for  a,  subject  to  a  <  1  and  then  for  a  using  the  following 
expression: 


bB  -  23c 
2aa  +  b 


-a  3)/a 


where 

•f  41  ->■  -V  -V 

a  =  2P  -P  -  2P  -P  -  2P  -P 
2  2  12  1  2 

B  =  2(P^  -P2  -P2  ) 


Note  that  in  deriving  the  above  expressions,  it  is  assumed  that  R1  is 
constant  but  R2  ,  in  general,  is  not. 

From  (1)  obtain 

xP  =  xi  +  a(x2-x  x) 

yp  =  y!  +  c<y2-yj) 

zp  =  Zi  +  a(z2-Zj) 

xp=  Xj  +  (x2-x  x)a  +  a(*,  -*2) 

yp=  ^i  +  (y2-yi  )a  +  <*(*2-y  x) 

zp=  Jx  +  (Zj-Zi  )&  +a(22-ij) 

Next  determine  the  latitude  and  longitude  of  the  point  P  and  rate  of 
change  of  these  quantities. 

<t>p  -  \  -  cos'1  (zp/Rp) 


Xp  =  tan~l(yp/xp) 


r 


♦  p  =  ip/  vx0  +'yfl 

•  _  xpj'p-yp^p 
X  p  =  Xp  +  y^> 

From  this,  determine  the  speed  and  direction  of  the  point  P  relative  to  a 
stationary  layer.  (Direction  is  relative  to  the  computation  system,  i.e., 
the  local  meridian  through  P.) 


%  -  cosOp  ) 

cos4>„ 


v  =  RpXp 
P  sinoip 

Next  determine  the  velocity  vx  (speed  and  direction)  of  the  point,  P, 
relative  to  the  moving  scintillation  layer.  This  is  given  by  the  vector 
vp  -  vs  where  vP  =  (vp,ap)  ,  vs  =  (vs,qs). 


vi  =  (vi  *°t  j )  where 

v,  =  [v£  +  v|  -  2vpvs  cos(otp-c(S)]1/2 


,  ,v0cosan-VcCOSOK 
a,=  tan-1  (-S-. — E — s.  ) 

'VnSiran-VflSinOs' 


Finally,  determine  the  velocity  of  P  in  the  local  coordinate  system  (see 
Figure  2). 
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APPENDIX  R  --  Computation  of  Theoretical  Spectra 

A  procedure  to  compute  the  theoretical  spectral  density 
function,  S(u))  (actually  S(f)  =  S(y£))  as  qi  ven  in  Section  V  is 
shown  below.  The  bulk  of  the  computations  is  done  by  device 
DISS,  PML.355  which  in  turn  uses  I M S L  subroutine  DCADRE,  to  per¬ 
form  the  required  integration.  The  parameters  used  in  this 
procedure  are  partitioned  into  two  classes  --  geometrical  and 
layer.  The  geometrical  parameters  are  those  which  define  the 
locations  and  motions  of  the  transmitter  (satellite)  and  receiver 
(aircraft).  (The  transmitter  frequency  is  also  included  among 
the  geometrical  parameters.)  Using  the  geometrical  parameters, 
the  necessary  "secondary"  parameters  can  be  computed  using  the 
algorithm  in  Appendix  A.  The  second  set  of  parameters,  the  layer 
parameters,  describe  the  turbulent  medium  in  terms  of  various 
correlation  function  coefficients  as  well  as  layer  height,  thick¬ 
ness  and  motion.  Refer  to  figures  3a,b,c  of  Section  V  and  figure 
2  of  Appendix  A  for  the  meaning  of  parameters  THETA,  PHI,  A,  PI 
and  P2.  Note  that  the  procedure  is  structured  so  that  the  user 
can  produce  several  spectra  at  one  run  by  varying  up  to  3  para¬ 
meters.  This  "family"  of  spectra  can  then  conveniently  be  dis¬ 
played  on  one  coordinate  frame  using  the  procedure  in  the  PLOT 
mode.  The  plots  are  produced  in  the  sane  manner  as  is  done  for 
the  experimental  spectra  (in  particular,  see  TM30  (Rarrett 
(1978))  for  an  example)  so  that  theoretical  and  experimental 
spectra  can  be  easily  compared. 
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APPENDIX  B  -  (Attachnent  1  ) 


A  Procedure  to  Conpute  Theoretical  Spectra 


.  PROC  ,  SCI  NTTH  , YEAR=7  6  ,  MONTH= 10  , DAY  =  20  ,  HOUP.=  1 ,  MI NUTE=8  , SECOND=0  , 

HA  =  3  0 , LAMBDA  =  $ -77  $,PHIA=$-11$,VA=470,ALPHAA=45,F=249, MOD EL =  2 , HS  =  300 , 
VS=70,ALPHAS=90 ,THETA=0 ,PHI=0,A=10,D=10,P1=3,P2=0, 

NN1=5, NN 2=1, NN 3=1, NP 1=11, NP 2=1, NP 3=1, 

P12=150,P13=$-81$,P14=90,P15=$-135$, 

P22=0,P23=0,P24=0,P25=0,P32=0,P33=0,P34=0,P35=0, 

PLOT=$-l$/$l . $,LIB=$JL/IN/SI/IM$ , I D= , I NPUTF=  #DATA , 

DEBUG  =  F/T ,SUP=0/1 .SWITCH l=OFF/ON , RED  I NK=OFF/ON . 

* 

.*  REVISION  —  JANUARY  29,1980 
.*  AUTHOR  —  BARRETT, TB 

.*  PURPOSE  —  RUN  (AND  DISPLAY)  THE  SPECTRAL  DENSITY  FOR  VARIOUS 
. *  THEORETICAL  LAYER  MODELS  AND  EXPERIMENTAL  GEOMETRIES. 

* 

. *  PARAMETERS  -- 

. *  (DEFAULTS  SHOWN  IN  PARENS,  THOSE  WITH  .  CAN  BE  GIVEN  AS  FLOAT) 


* 

(1) 

YEAR 

(76) 

YEAR  OF  EXPERIMENT 

it 

(2) 

MONTH 

(10) 

MONTH  OF  " 

★ 

(3) 

DAY 

(20) 

DAY  OF 

★ 

(4) 

HOUR 

(1) 

HOUR  OF  " 

★ 

(5) 

MINUTE 

(8) 

MINUTE  OF  " 

★ 

(6) 

SECOND 

(0.  ) 

SECOND  OF  " 

★ 

THE  ABOVE  PARAMETERS 

ARE  USED  TO  OBTAIN  SATELLITE  DATA 

★ 

(7) 

HA 

(30.  ) 

AIRCRAFT  HEIGHT  (THOUSANDS  OF  FEET) 

★ 

(8) 

LAMBDA 

(-77. ) 

AIRCRAFT  LATITUDE  (DEGS),  -  IS  SOUTH  OF  EQ. 

★ 

(9) 

PH  I A 

(-11. ) 

AIRCRAFT  LONGITUDE  (DEGS),  -  IS  WEST  OF  GR. 

★ 

(10) 

VA 

(470. ) 

AIRCRAFT  GROUND  SPEED  (KNOTS) 

* 

(11) 

ALPHAA 

(45.  ) 

AIRCRAFT  HEADING  (DEGS  CLOCKWISE  OF  NORTH) 

★ 

(12) 

F 

(249. ) 

TRANSMITTER  FREQUENCY  (MHZ) 

* 

LAYER 

PARAMETERS  — 

★ 

(13) 

MODEL 

(2) 

1  =  >GAUSS I AN , 2  =  >POWER  LAW , 3  =  >KOLMOGOROF 

* 

(14) 

HS 

(300  .  ) 

LAYER  HEIGHT  (KM) 

★ 

(15) 

VS 

(70.  ) 

LAYER  SPEED  (METERS/SEC) 

* 

(16) 

ALPHAS 

(90.  ) 

LAYER  DIRECTION  OF  MOTION  (DEG  WRT  GEO.  N.) 

★ 

(17) 

THETA 

(0.  ) 

VERTICAL  TILT  OF  "CORRELATION  BLOB"  (DEGS) 

★ 

(18) 

PHI 

(0.  ) 

"SLEW  ANGLE"  (DEGS)  OF  THE  "CORRELATION  BLOB" 

★ 

FROM  THE  "LOCAL"  X-AXIS. 

★ 

(19) 

A 

(10.  ) 

CORRELATION  BLOB  AXIS  RATIO  (>1) 

it 

(20) 

D 

(10.  ) 

LAYER  THICKNESS  (KM) 

it 

(21) 

PI 

(4.  ) 

CHARACTERISTIC  CORRELATION  LENGTH  OR  OUTER 

* 

SCALE  LENGTH  (KM) 

★ 

(22) 

P2 

(1.  ) 

INNER  SCALE  LENGTH  (METERS) 

* 

THE  FOLLOWING 

PARAMETERS  PERMIT  VARIATIONS  IN  UP  TO  3  PARAMETERS. 

* 

NN1 

(5) 

TOTAL  NUMBER  OF  1ST  PARAMETER  (1  TO  5) 

* 

NN  2 

(1) 

TOTAL  NUMBER  OF  2ND  PARAMETER  (1  TO  5) 

★ 

NN  3 

(1) 

TOTAL  NUMBER  OF  3RD  PARAMETER  (1  TO  5) 
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* 

NP1 

(ID 

PARAMETER  TO  BE  VARIED 

* 

NP2 

(1) 

ETC.  FOR  2ND  AND  3RD  PARAMETERS. 

* 

NP3 

(1) 

* 

PI  2 

(150.  ) 

VALUE  OF  FIRST  VARAI ATION  OF  FIRST  PARAM. 

* 

Pl  3 

(-81.) 

SECOND 

* 

PI  4 

(90.  ) 

THIRD 

* 

PI  5 

(-135.) 

FOURTH 

* 

P22  - 

P25 

* 

P32  - 

P35 

SIMILARLY  FOR  THE  OTHER  VARIED  PARAMETERS. 

* 

THESE  ARE  ALL  DEFAULTED  TO  0. 

* 

LIB 

( $JL/IN/SI/IM$ )  LIBRARY  LOCAL  FILE  NAMES  TO  BE  USED. 

* 

PLOT 

(-1/$1.$ 

)  IF  PLOT  IS  SPECIFIED  PLOTS  ARE  PRODUCED. 

* 

ID 

(NULL) 

IF  SPECIFIED,  TAPE3  FROM  DISS  IS  CATALOGED 

* 

UNDER  THIS  ID  (AS  TAPE3,ID=ID) 

* 

DEBUG 

(F/T) 

IF  SPECIFIED,  DEBUG  INFO  IS  PRODUCED. 

* 

SUP 

(0/1) 

IF  SPECIFIED,  IMSL  ERRORS  GOTO  TAPEl 

* 

SWITCH 1  (OFF/ON) 

IF  SPECIFIED,  PLOTS  ARE  SINGLE  COLOR. 

★ 

REDINK 

(OFF/ON) 

IF  SPECIFIED,  PLOTS  ARE  IN  RED  INK 

( SWITCH 1  SHOULD  ALSO  BE  ON) 


*  ******************************************************************** 
* 

*  ***************************£  x  A  M  P  L  E  ************************* 

* 

*  BARTH, CM120000 ,T511 . 

*  ATTACH (JL,JURGLIBX3693818,ID=BANDES) 

*  ATTACH (IN, INFOLIBX3693818,ID=BANDES) 

*  ATTACH (SI ,SIGLIBX 3693818, ID=WONG) 

*  ATTACH (IM, I MS LLIBX 3693818, ID=WONG ) 

*  LIBRARY (JL) 

*  ATTACH (SATFIL,Y76M10D20X3693818 ,ID=BANDES) 

*  SCI NTTH,ID=BARRETT, PLOT. 

* 

*  ******************************************************************* 


* 


IFE , $SWITCHl$=$ON$ , LAB0 . 

SWITCH (1) 

ENDIF,LAB0 . 

REQUEST (TAPE3,*PF) 

FTN (I=INPUTF,SYSEDIT,R=3,SL,L=Q) 
LDSET (#LIB=LIB) 

LGO ( PL=60000 ) 

EXIT (U) 

IFE , $DEBUG$=$T$ , LAB3 . 

REVERT. 

ENDIF, LAB3 . 

RETURN (LGO) 

IFE, $ID$.NE.$$, LABI. 

CATALOG (TAPE 3,# ID- ID) 
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ENDIF, LABI . 

I FE , $  PLOT? . EQ . $ 1 . $ , LAB2 . 

ATTACH (PEN,ONLINEPEN) 

ATTACH ( PLT , PLOTL IB ) 

LIBRARY (PEN, PLT) 

REQUEST ( #PLOT , *Q) 

IFE,$REDINK$=$OFF$,LAB4. 

DISPOSE ( #PLOT, *OL) 

ELSE , LAB4 . 

DISPOSE (#PLOT, *PL ) 

ENDIF, LAB4. 

ENDIF, LAB2. 

FTN ( I=INPUTF,SYSEDIT,R=3,SL ,L=Q) 

LDSET ( PRESET=0 , #LIB=LIB ,SUBST=GENGENU-GENGEN4) 

LGO (PL=60000 ) 

REWIND (TAPE2) 

COPY (TAPE2 ) 

REWIND (Q) 

COPY (Q) 

EXIT. 

DMP ( 60000 ) 

REWIND (TAPE2) 

COP  Y ( TAPE2 ) 

REWIND (Q) 

COPY ( Q ) 

.DATA 

PROGRAM  RUNDIS (OUTPUT  =  6  4 , TAPE 3 ,SATFIL,TAPEl=64) 
COMMON/GEOLAY/GL (22,5) 

COMMON/PH I AUX1/# DEBUG 
LOGICAL  #DEBUG 
DIMENSION  IGL (22,5) 

EQUIVALENCE  (GL, IGL) 

#DEBUG=. DEBUG. 

NY=NNl*NN2*NN3 
DO  90  1=1, NY 
IGL ( 1 , I ) =YEAR 
IGL (2, I) =MONTH 
IGL (3,1) =DAY 
IGL (4,1) =HOUR 
IGL (5,1) =MI NUTE 
GL (6 , I ) ^SECOND 
GL ( 7 , I ) =HA 
G  L ( 8 , 1 ) = LAMBDA 
GL(9, I) =PH I A 
GL(10, I) = VA 
GL ( 1 1 , I) =ALPHAA 
GL ( 1 2 , I ) =F 
IGL (13,1) =MODEL 
GL ( 1 4 , I ) =HS 
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GL (15,1) =VS 
GL (16 , I ) =ALPHAS 
GL (17,1) =THETA 
GL (18,1) =PHI 
GL (19,1) = A 
GL (20,1) =D 
GL (21,1) =Pl 
GL (22,1) =P2 
CONTINUE 

IF  (NN 1  .EQ.  1)  GOTO  600 

IF  (  (NPl  .LT.  6)  .0.  (NP1  .EQ.  13))  GOTO  100 

GL(NP1 ,2) =P12 

GL (NPl , 3 ) =P1 3 

GL (NPl , 4 ) =P1 4 

GL(NP1 ,5) =P15 

GOTO  200 

IGL (NPl , 2 ) =P1 2 

IGL (NPl , 3 ) =P1 3 

IGL (NPl , 4 ) =P1 4 

IGL (NPl , 5 ) -Pi  5 

IF  (NN2  .EQ.  1)  GOTO  600 

IF  (  (NP2  .LT.  6)  .0.  (NP2  .EQ.  13))  GOTO  300 

GL (NP2 , 2 ) =P2  2 

GL (NP2 , 3 ) =P2  3 

GL (NP2 , 4 ) =P2  4 

GL(NP2 ,5) =P25 

GOTO  400 

IGL (NP2 , 2 ) =P22 

IGL (NP2 , 3 ) =P2  3 

IGL (NP2 , 4 ) =P2  4 

IGL (NP2 , 5) =P25 

IF  (NN 3  .EQ.  1)  GOTO  600 

IF  (  ( NP 3  .LT.  6)  .0.  (NP3  .EQ.  13))  GOTO  500 

GL (NP3 , 2) =P32 

GL (NP3 , 3 ) =P33 

GL (NP3 , 4 ) =P34 

GL (NP3 , 5 ) =P3  5 

GOTO  600 

IGL (NP3 , 2)  =P3  2 

IGL (NP3 , 3 ) =P33 

IGL (NP3 , 4 ) =P3  4 

IGL (NP3 , 5 ) =P3  5 

CALL  D ISS ( NN 1 , NN 2  ,  NN  3  , NP 1 , NP 2 , NP 3 , SUP ) 

END 

PROGRAM  TESTGEN (OUTPUT=64 ,TAPE2=64 ,TAPE3) 
COMMON  DATA (7000) 

COMMON/CONTROL/IC (100) 

COMMON/ AUX I L/I A (10) 


COMMON/PROCON/IPR ( 700 ) 

COMMON/PRI NT/IPRNT ( 300 ) 

COMMON/FIT/FIT (600) 

COMMON/SUB/SUB , Ml , Ll , I WHER ( 5 ) 

COMMON/M I NPLT/M I NPLT (10) 

COMMON/GENGENA/I WHERE 
DIMENSION  PRNT ( 300 ) 

EQUIVALENCE  ( IPRNT, PRNT) 

DATA  (IPR=256, 100. ,9. ,2,6,50) 

DATA  { I PRNT=PLOT ,0,1, "LOG" ,-1 . ,4. ,0., 1,0,1, "LOG" ,-8. ,1. ,1 

IWHERE=3000 

NY=NNl*NN2*NN3 

IPRNT (2) = I PRNT (9)  =NY 

IF  (NY  .EQ.  1)  GOTO  100 

K=  15 

DO  110  1=2, NY 

IPRNT (K) =1 

PRNT (KFl) =PRNT (11) 

PRNT (K+2) =PRNT (12) 

PRNT ( K+3 ) =PRNT (13) 

PRNT ( K+4 ) =PRNT (14) 

K  =  K  +  5 

110  CONTINUE 

100  CALL  GENGEN (1 , 1 ,1 ,270,1 ,0,NVAL) 

IF  ( IPR ( 7  )  .EQ.  0)  GOTO  100 

CALL  FINISH 

END 
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APPENDIX  C  --  Computation  of  the  Spectral  Density  Function 
Estimates 


The  estimator  which  was  used  is  that  given  by  equation  7.7 
in  Section  VII.  For  technical  reasons  we  write  this  equation  in 
the  form: 

1  2mf|  -1 

(1)  SN(w)  =  j- -  1  w(n/nN)RN(n)e-ino. 

(Note  that  2m^  is  the  actual  truncation  point  rather  than  m^ 

because  w  has  been  "normalized"  to  be  0  for  |  x  |  2l2  rather  than 
for  |x|  >1  as  in  Equation  (7.7).  More  will  be  said  about  w 
later.)  Let  C|^(n)  =  w(n/m^)  •  R^(n)  and  consider  S^  evaluated  at 
the  points 

^  j  /  ^  |\|  j-D»***.2mj(j-l. 


Then : 

(2) 


1  2mN-l 

J: _  £ 

2  n  =  - 2m^  +  1 


CN(n)e-i^nJ'/mN 


I2en_1 

*n  =  0 


CN(n)cos(^) 


CN(0) 
2  * 
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Note  that  since  is  real.  Re  {OFT^C^  }  is  even  (see  Appendix  A 
of  Barrett  (1973))  and  hence  needs  to  be  evaluated  only  at  the 
points  j  =0  , .  .  .  ,  my/ 2  =  m^  .  Of  course,  using  the  FFT  algorithm  S^j 
is  found  at  all  points  j  =0  , . . . ,  mj- 1  and  Im  {DFTm^.C^  }  is  computed 
whether  one  wants  it  or  not.  We  are  concerned  mainly  with  evalu¬ 
ating  at  the  points  nj/m^  because  of  correlation  considerations 
mentioned  in  Section  VII.  It  is  interesting  to  note,  however, 
that  knowledge  of  R e  tDFTm^C^  }  and  Irn  {DFT^^Cfj  }  is  necessary  and 
sufficient  for  evaluating  over  the  continuum  f  0 ,  it "!  because  of 
the  "sampling  theorem"  (see  Appendix  A  of  Barrett  (1973)). 

In  this  next  section  we  will  be  concerned  mainly  with  the 
computation  of  the  correlograms  R^.  First,  remember  that  the 
theoretical  properties  of  this  class  of  estimators  are  based  on 
their  asymptotic  behavior  as  N,  the  length  of  the  time  series  (or 
"amount  of  data")  N  becomes  larger  and  larger.  Second,  note  that 
the  truncation  point  mj  must  grow  larger  at  a  much  smaller  rate 
than  N  for  good  asymptotic  behavior  of  S^.  This  means  that  it  is 
not  necessary  to  compute  the  correlogram  R  ^ ( n )  for  all  possible 
values  of  n,  and,  in  fact,  as  N  becomes  large,  need  be  cal¬ 
culated  over  only  a  relatively  small  percent  of  its  support.  As 
will  be  shown,  it  is  convenient  to  calculate  over  a  fixed 
interval  (say  0,...,Nj-l),  where  is  equal  to  or  greater  than 
the  largest  value  of  the  truncation  point,  regardless  of  the 
amount  of  data  available.  (It  goes  without  saying  that  in  any 
finite  experiment,  the  amount  of  data  and  hence  the  truncation 
point  will  presumably  be  the  truncation  point  used  in  the  last 
"trial"  (see  below).)  The  reason  for  this  is  that,  1)  the  FTT  is 
a  very  efficient  algorithm  for  computing  correlations  if  is 
"highly  composite"  (in  the  algorithm  used  this  means  Ni  =  2'<>  k  a 
positive  integer),  and  2)  it  is  not  necessary  to  recompute  "all 
of"  R^j  as  more  data  is  added.  In  fact,  consider  a  time  series  of 
length  N  such  that  N=£N^  where  l  is  a  positive  integer.  Then  it 
is  apparent  (see  Figure  Cl.)  that  R^(n')  can  be  written  as: 

RN  (  n  '  )  =  ^ik|1rrk(n')  +  crk  (  n  ' )  J  =  R  (n1)  , 

n ' =0  , . .  .  ,Nj-l 


66 


0  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15 


0  1  2  3  4  5  6  7 
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*  ->  R 16  ( 3 ) 


->  rrk( 3) 


->  rrk  +  i ( 3 ) 


*  ->  crk (3) 


*  means  multiply  corresponding  data  then  add.  Corresponding  data  are  shown 


as  .  in  the  same  columns. 


R 16 ( 3 )  E  R2  ( Nx  )  <  3  )  =  [rrk(3)  +  rrk+i(3)  +  crk(3)] 

Figure  Cl.  "Decomposition"  of  the  correlogram  for  Nk=8, 
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where 


N  !  -  n  '  - 1 

rrk(n')  =  Ln  xk  (n)xk  (n  +  n  '  ) 


n  '  -1 

crk(n')  =  Ln  xk(n+N1-n'  )xk  +  i(n)  , 
n  =  u 

xk  is  the  kth  block  of  data  and  the  data  in  each  block  is 
numbered  from  0  to  Nj-1.  Moreover,  rrk  and  crk  can  be  found  by 
"circular"  correlation  and  thus  by  a  double  discrete  Fourier 
transform  as  will  be  described  in  some  detail  below.  In  fact  if 

\ 

f(n)  =  xk(n)  n  =  0  , . .  .  ,Nj- 1 

=  0  n-Nj . 2  N  !  - 1 


g( n)  =  xk+1 (n) 


n=*0 , . . .  ,Nj - 1 


=  0 


naNj,...,2Nj-l 


then  rrk(n‘)  =  f*f(n')  n  '  =  0  , .  . .  ,N j - 1 

crk(n')  =  f*g(Nj+n‘)  n ' =0  , . .  .  ,  N j - 1 

where  *  means  circular  correlation  for  G  =  Z i  (see  Appendix  A  of 
B  a  r ret t  (  1  973  )). 

In  obtaining  the  spectral  density  function  estimates,  two 
"devices"  are  used.  (A  description  of  these  "devices"  in  the 
form  of  comment  lines  from  the  source  software  is  given  as 
attachments.)  The  first  one  is  called  CORRFL  and  computes  the 
required  correlogram  R^(n)  from  scintillation  data  input.  The 
second  device,  called  SPECTR,  uses  the  output  of  CORREL  to 


68 


compute  the  spectral  density  function  estimator,  S^.  For 
technical  reasons  (storage  considerations  and  run  time)  it  is 
convenient  to  make  several  passes  (these  were  referred  to  as 
"trials",  l ,  above)  through  CORREL,  obtaining  additional  "blocks" 
of  scintillation  data  each  time  until  R^  has  been  computed  for 
all  required  lags.  Following  this,  SPECTR  is  used  to  calculate 
S^.  An  example  of  this  process  can  be  found  in  Barrett  (  1978  ), 
p.  27-33.  CORREL,  in  particular,  utilizes  many  of  the  symmetry 
properties  of  the  FFT. 

The  computat i ons  are  outlined  below.  Figure  C2  is  an  array 
diagram  summarizing  the  major  computat i onsl  subdivisions.  Assume 
that  the  (A-l)th  correlogram  has  been  found  and  that  the  discrete 
Fourier  transform  Xj£_i  of  the  corresponding  block  of  data  is 
stored.  (In  the  sequel,  refers  to  the  discrete  Fourier 
transform  of  the  At h  block  of  data  o_r  this  transform  in  reverse 
binary  order.  When  computing  correl ograms ,  it  is  not  necessary 
to  put  transforms  in  correct  order;  the  FFT  algorithm  used 
naturally  outputs  in  reverse  binary  order.)  Then  two  new  blocks 
of  data  and  x^  +  i  are  produced  as  described  in  Barrett  (1978) 
to  be  used  as  "input"  to  CORREL.  Zeroes  are  appended  to  these 
data  blocks  to  obtain  the  equivalent  of  functions  f  and  g 
described  above.  (The  same  symbol  x^  will  denote  x^  or  x^ 
appended  with  zeroes;  similarly  X ^  will  denote  the  Fourier 
transform  of  x^  appended  with  zeroes.)  The  two  real  "functions" 
xA  and  X£+i  are  simultaneously  transformed  as  described  in 
Appendix  A  (Barrett  (1973))  -  (see  Fact  A).  Then  X^X^and  X^+jX^+i 
are  formed  and  inverse  transformed,  again  simultaneously,  to  ob¬ 
tain  the  quantities  rr  A  and  rr^+^.  rr4  is  added  to  Rjj,-1  •  ty  ( *-l  ) 
and  rr^+j  is  stored.  Next  the  quantity  X^.^T^  is  formed  and 
inverse  transformed  to  obtain  crA_i  and  RA=  R*-i  V ( *-1 )+rr^+cr^_i 
output.  (In  performing  the  inverse  transfrom.  Fact  B  of  Appendix 
A  (Barrett  (1973))  is  used.)  Finally  X^X^+i  is  formed  and 
inverse  transformed  and  R^+i  found  from  R*+i  =  RA*V(*  +  1)  +  rr*+i 
+  cr^  to  complete  one  cycle  through  CORREL. 
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The  correlogram  computations  are  based  on  Comm.  A.C.M. 
algorithm  345,  an  Algol  "convolution"  procedure  by  Singleton 
(1969a,  1969b).  This  algorithm  consist  essentially  of  four 
procedures:  FFT4  produces  the  discrete  Fourier  transform  of  a 
complex  valued  function  (in  reverse  binary  order);  REVFFT4 
produces  the  inverse  DFT  of  complex  data  which  is  in  reverse 
binary  order;  REALTRAN  unscrambles  the  transform  (or  inverse 
transform)  if  Fact  C  (Appendix  A  of  Rarrett  (1973))  is  used  to 
find  the  transform  of  a  real  function  of  even  length; 
REVERSEBINARY  (which  is  not  used)  is  used  to  return  the  output  of 
FFT  to  normal  order.  Using  this  nomenclature,  CORREL  can  be 
summarized  as  follows: 

1)  x  and  xfc+i  are  the  inputs  to  FFT4  which  outputs  X  ^  and  XA+i 
in  scrambled  form.  This  scrambled  output  is  unscrambled  as 
described  in  Fact  A  to  obtain  X  ^  and  X^+j  which  are  stored. 

2)  X  Y£  and  XA  +  1'XJl+i  are  used  as  inputs  to  REVFFT4  which  returns 
the  scrambled  versions  of  rr^  and  rr^+j.  These  are 
unscrambled  as  above. 

3)  X^iXjj  are  formed  and  then  the  real  function  Re  [  X  £_  i  X  z  ]- 
ImCX^.jX^].  This  latter  function  is  nade  into  the  real  and 
imaginary  parts  of  a  complex  function  which  in  turn  is  inverse 
transformed  by  REVFFT4  and  "unscrambled"  next  by  REALTRAN  to 
obtain  cr j . 

4)  X^X^+i  is  formed  and  cr  found  as  above. 

Following  CORREL,  the  spectral  density  function  estimate  is 
found  by  SPECTR.  First  a  subroutine  (which  nay  be  changed  for 
different  lag  windows)  produces  the  weighted  correl  ograms ,  i.e., 
the  Cfj(n)  referred  to  in  Equation  (2).  The  discrete  Fourier 
transform  of  CN  is  performed  using  a  fast  Fourier  transform 
algorithm  provided  by  Norman  Brenner  formerly  of  MIT  Lincoln 
Laboratory.  This  algorithm  may  be  used  with  data  of  any  length 
although  it  works  faster  with  a  highly  composite  length.  The 
original  program  may  be  used  to  transform  real  or  complex  valued 
functions  defined  on  x  ZN?  x  ....  ZNn  where  Nj ,N2 , . . • ,Nn  may 
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be  any  integer.  SPECTR  uses  an  abridged  version  of  this 
algorithm  which  transforms  real  functions  on  where  Nj  is 
even . 

The  lag  window  w  which  was  used,  has  some  nice  properties 
and,  in  fact,  was  chosen  because  of  these  properties.  (See 
Woodroofe  and  VanNess  (1967).)  The  Fourier  transform  of  w,  the 
"spectral  window"  W,  is  assumed  to  be  non-negative  even  and 
continuous  with  a  continuous  Fourier  transform  W  =  w  (the  "lag 
window").  In  addition,  we  assume  that  w  has  the  following 
propert i es : 


(a) 

w  ( x  )  =  0 

|x  |  >2 

(b) 

w( 0 )  =  1 

(c) 

W  *  (0) 

exists. 

The  assumptions  for  W  and  assumption  (b)  and  (c)  for  w  are 
conditions  imposed  by  theorem  2.2  of  Woodroofe  and  VanNess 
(1967).  Condition  (a)  is  imposed  partly  for  computational 
convenience  and  partly  because  theorem  2.2  is  based  on  the  fact 
that  truncated  spectograph  estimators  are  assumed.  Normally,  w 
is  truncated  at  |x  |  =  1;  for  technical  reasons  truncation  at  lx  I 
=  2  is  used. 

A  function,  w,  which  satisfies  these  conditions  is  that 
obtained  by  self-convolving  a  suitable  "triangle"  function. 

This  particular  lag  window  was  apparently  first  proposed  by 
E.  Parzen  in  1957  (see  Parzen  (1961)).  It  is  probably  the 
simplest  example  of  a  function  which  meets  the  conditions. 

The  Fourier  transform  of  such  a  function  is  obviously 
non-negative  since  the  Fourier  transform  of  the  triangle  function 
is  non-negative.  Let  the  triangle  function  be  given  by: 

Mt)  =  a(l  -  \^\  It  |  <  t 

=  0  |t  |  >  T  . 
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Then  the  function  obtained  from  sel f -convol ut i on  of  this  function 
i  s : 

A*A(t)  =  a2i2(\(±)3)  -  (4)2  +  |)  0  <  t  <  tu 

=  (2  -  y)3  T  <  t  <  2  T 

=  n  t  >  2  T 


and  by  even  symmetry  for  t  <  0. 

If  t  =  1  and  ,.2  =  3/2,  then  we  obtain 

w ( t )  =  |t3  -  |t2  +1  0  <  t  <  1 

=  |(2  -  t)3  1  <_  t  <  2 

=  0  t  >  2  . 
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X1 
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A  of  Barrett  ( 1973) . ) 

REVFFT4  plus  Fact  A 

X  1  X 
CO  j  CO 

1 

1 

! 

1 

Q 
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X1 
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Q  =  Rj  +  rr2 
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x3 
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Ro  =  yRi  +  rr2  +  cr^ 
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S3  output 

Fact  B  and  Fact  C 

x2*3  X2X3 
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(1)  even  numbered  X2X3 

(1)  (2) 

rr3 

x3 

(2)  odd  numbered  X2X3 

REVFFT4  and  REALTRAN 
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2 

R3  =  ^2  +  rr3  +  cr2 

If  SPECTR  called  here 

s3 

r3 

x3 

x3 

(9) 

S3  output 

(Return  to  (1)) 

Figure  C2. 

—  Array 

use 

diagram  of  CORREL  and  SPECTR. 

(Note  that 
some  stages 

each  array 
each  half 

is  shown  in  t wo  halves  since  in 
is  used  to  store  different  data 

73 


APPENDIX  C  -  (Attachment  1) 


Comments  from  CORREL,  the  Correloqram  Estimator 


SUBROUTINE  CORREL ( L , I A , I X , I Y , I PRT , I PLT , I NDIC , NSI ZE , 
$  A,B,C,D,Q1 ,  Q3 , S I NE , IGO) 


C 

c 

c 

c 

c 
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PML  — 
REVISION  — 
AUTHOR  -- 
PURPOSE  — 


236 

MARCH  6,1979 
BARRETT ,TB 

COMPUTES  THE  CORRELOGRAM  FOR  A  REAL  TIME  SERIES 
WHCIH  IS  GENERATED  IN  BLOCKS  OF  SIZE  IBLOCK  SUCH  THAT 
I  BLOCK  IS  A  POWER  OF  2  AND  ( IBLOCK-1 ) *IA  (1 )  IS  THE 
MAX.  LAG  VALUE.  THIS  DEVICE  IS  SIGPRODS 
COMPATIBLE  WITH  STATE  STORED  EXTERNALLY,  LARGELY 
IN  ARRAY  A  ->  SINE. 


(PARTIALLY) 


USE  OF  STANDARD  ARGUMENTS 

L  START  OF  /CONTROL/  TO  STORE  •’STATE" 

AS  FOLLOWS 

1  INITIAL  CALL  INDICATOR 

2  NUMBER  OF  OUTPUT  TERMINALS  ON  INPUTTING 

DEVICE 

3  STARTING  LOCATION  IN  DATA  FOR  OBTAINING  INPUT 

4  STARTING  LOCATION  IN  DATA  FOR  STORING  OUTPUT 

5  I I =NO .  OF  BLOCKS  WHICH  HAVE  BEEN  PROCESSED 

6  MZ, ( 2 * *M Z ) =NZ = 2* I BLOCK 

IA  START  OF  /AUXIL/  TO  OBTAIN  AUXILIARY  INPUT  DATA 
AS  FOLLOWS 

1  TIME  INTERVAL  TO  BE  ASSOCIATED  WITH  1  LAG  UNIT 

2  UNITS  (HOLLERITH)  TO  BE  ASSOCIATED  WITH  THE  TIME  INTERVAL 
THIS  IS  STORED  AS  PART  OF  THE  PIN  LABEL  FOR 

LABELING  PURPOSES  (ONLY  USED  IF  INDIC=2) 

IX  INPUT  TERMINAL  LOCATION  (NOTE  THIS  IS  A  SINGLE 
INPUT  (REAL)  DEVICE) 

IY  OUTPUT  TERMINAL  (CAN  BE  AN  ARRAY  OF  LENGTH  2) 

CONTAINING  THE  LOCATIONS  FOR  UP  TO  2  OUTPUTS. 

THE  FIRST  IS  THE  MAIN  OUTPUT  (CORRELOGRAM) 

WHILE  THE  SECOND  IS  THE  "TIME  BASE"  OUTPUT  IF 
I ND  IC=  2 

IPRT  PRINTER  NUMBER 

IPLT  RECORDER  NUMBER 

USE  OF  PARTICULAR  ARGUMENTS 

INDIC  (1,2, 3, 4)  INDICATING  THE  MODE  OF  OPERATION 
IF  I ND IC= 1  ONLY  THE  CORRELOGRAM  IS  PRODUCED. 

IF  I ND IC  =  2  A  TIME  BASE  OUTPUT  AND  CORRELOGRAM  ARE  PRODUCED 
IF  I NDIC=  3  A  TIME  BASE  OUTPUT  AND  NORMALIZED  (TO  1) 

CORRELOGRAM  IS  OUTPUT 

IF  I ND IC=  4  A  NORMALIZED  CORRELOGRAM  AND  NO  TIME  BASE 
IS  PRODUCED. 

NSI ZE  IS  THE  SIZES  OF  THE  ARRAYS  A->Q3. 
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NSI  ZE 


NSI  ZE 


C  SHOULD  BE  AT  LEAST  2*MT+2 

C  A- >Q3  ARRAYS  FOR  STORING  MOST  OF  THE  STATE  OF  CORREL 

C  SINE  ARRAY  FOR  STORING  VALUES  OF  FUNCTION  SIN 

C  (MIN.  SIZE  M2  FROM  2*MT=2**MZ) 

C  IGO  SEE  OTHER  OUTPUT 
C  OTHER  INPUT  - 

C  /PROCON/1  WHICH  IS  THE  DATA  BLOCK  SIZE  IS 

C  ASSUMED  TO  BE  THE  VALUE  OP  MT 

C  /PROCON/8  IS  THE  ACTUAL  AMOUNT  OF  DATA  GENERATED 

C  (NOTE  THAT  THE  LAST  BLOCK  MAY  NOT  BE  COMPLETELY 

C  FILLED  WITH  DATA) 

C  OTHER  OUTPUT  - 

C  IGO  (INDICATOR  OUTPUT  FOR  USE  BY  OTHER  DEVICES  AS  FOLLOWS) 

C  0  =>  NO  CORRELOGRAM,  RETURN  TO  GENERATOR 

C  1  =>  CORRELOGRAM,  RETURN  TO  CORREL 

C  2  =  >  CORRELOGRAM,  RETURN  TO  GENERATOR 

C  (RETURNS  TO  GENERATOR  ARE  CONDITIONAL  ON  IPP ( 7 ) =0 ) 

C  SUBROUTINES  (NON-SYSTEM) 

C  FFT4  FAST  FOURIER  TRANSFORM 

C  RVFFT4  INVERSE  FFT 

C  REALTR 

C  SPECIAL  PRECAUTIONS  - 

C  1  /PROCON/1  SHOULD  BE  A  POWER  OF  2.  IF  IT  IS  NOT 
C  THE  SUBROUTINE  WILL  TERMINATE.  NOTE  HOWEVER 

C  THAT  THE  DATA  BLOCK  DOES  NOT  HAVE  TO  BE  COMPLETELY 

C  FILLED 

C  2  THE  ARRAY  A->Q3  MUST  BE  DECLARED  IN  THE  CALLING 
C  PROGRAM  WITH  SUFFICIENT  SIZE  (SO  THAT  2  DATA 

C  BLOCKS+2 )  WILL  FIT  INTO  EACH  ARRAY. 
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APPFNDIX  C  -  (Attachment  2) 


Comments  from  SPECTR, 


The  SPfir  Estimator 
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SUBROUTINE  SPECTR (L , I A , IX, IY , IPRT , IPLT ,TL , IGO,A,D) 

PML  231 

REVISION  --  JUNE  2,  1977 
AUTHOR  —  UNKNOWN 

PURPOSE  —  COMPUTES  THE  SPECTROGRAM  OF  A  REAL  TIME  SERIES  FROM  THE 
CORRELOGRAM  AS  PRODUCED  BY  CORREL.  THIS  DEVICE  IS  SIGPRODS 
COMPATIBLE  WITH  STATE  STORED  EXTERNALLY.  IT  HAS  TWO  OUTPUT 
TERMINALS  THE  FREQUENCY  BASE  AND  THE  TIME  BASE.  NOTE  THAT 
THIS  IS  CONSIDERED  TO  BE  ESSENTIALLY  A  TERMINATING  DEVICE 
AND  IS  NOT  IN  THE  "MAIN"  SIGNAL  FLOW  STREAM. 

USE  OF  STANDARD  ARGUMENTS 

L  START  OF  /CONTROL/  TO  STORE  "STATE" 

1  INITIAL  CALL  INDICATOR 

2  ACTUAL  "TIME"  INTERVAL 

3  MT  FOR  PREVIOUS  CALL 

IA  START  OF  /AUXIL/  TO  OBTAIN  AUXILIARY  INPUT  DATA  AS  FOLLOWS 

1  "TIME"  INTERVAL  ASSOCIATED  WITH  THE  CORRELOGRAM 
(FLOATING  POINT).  IF  NOT  LEGIT-  “'ATE  IS  SET  TO  1. 

2  "TIME"  INTERVAL  UNITS  (HOLLERITH).  (NOTE  THAT 

THE  UNITS  SHOULD  BE  8  CHARACTERS  OR  LESS  -  LEFT 

JUSTIFIED 

THE  MAXIMUM  FREQUENCY  VALUE  GENERATED  BY  SPECTR  IS 
1 . / (2 . *"TIME"  ) . 

IX  INPUT  TERMINAL  LOCATION  IN  DATA  (INPUT  MUST  BE 
CORRELOGRAM  OUTPUT 

IY  OUTPUT  TERMINAL  (ARRAY  OF  LENGTH  2) 

THE  FIRST  ELEMENT  REFERS  TO  THE  LOCATION  IN  DATA  FOR 
THE  SPECTROGRAM  WHILE  THE  SECOND  LOCATES  THE  FREQUENCY 
BASE  OUTPUT. 

IPRT  PRINTER  NUMBER 

IPLT  RECORDER  NUMBER 

USE  OF  PARTICULAR  ARGUMENTS 

TL  TRUNCATION  TIME, MUST  BE  LESS  THAN  IBLOCK*TAU  WHERE  TAU 
IS  THE  SAMPLE  INTERVAL  (=AI(IA)  ) 

IGO  =0  DOES  NOTHING 

=1  COMPUTE 

THIS  INDICATOR  IS  NECESSARY  BECAUSE  OF  THE  WAY  THAT 
CORREL  INTERACTS  WITH  THE  DATA  SOURCE. 

A , D  WORK  STORAGE  ARRAYS  -  USUALLY  THE  SAME  ONES  USED  BY 
CORREL.  (SEE  CORREL) 

OTHER  INPUT  - 

/PROCON/1  WHICH  IS  THE  DATA  BLOCK  SIZE,  IBLOCK,  IS 

ALSO  USED  TO  DETERMINE  THE  AMOUNT  OF  DATA  OUTPUT  (AS 
WELL  AS  INPUT).  THE  NUMBER  OF  "FUNDAMENTAL"  AND  UNIQUE 
FREQUENCIES  IS  MT/2+1  WHICH  AT  MOST  IS  HALF  OF  IBLOCK. 

THE  OTHER  FREQUENCIES  ARE  FOUND  BY  INTERPOLATION  USING 
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THE  SAMPLING  THEOREM. 

OTHER  OUTPUT  - 

NONE 

OTHER  SUBROUTINES  (NON-SYSTEM) 

FFT  GENERAL  PURPOSE  FFT  ALGORITHM  FOR  REAL  DATA 
ALSO  DOES  INTERPOLATION) 

W8AUTO  WEIGHT  FUNCTION  FOR  WEIGHTING  THE  CORRELOGRAM 

SPECIAL  PRECAUTIONS 

1  SPECTR  CAN  ONLY  BE  ATTACHED  TO  CORREL  OR  POSSIBLY  A 
"RECORDER"  WHICH  HAS  STORED  THE  OUTPUT  OF  CORREL. 

2  THE  TRUNCATION  POINT  MUST  BE  LESS  THAN  OR  EQUAL  TO 
IBLOCK  (IPR(l)).  IF  IT  IS  NOT  THE  PROGRAM  WILL 
SET  IT  TO  THIS  VALUE. 
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